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I .   INTRODUCTION 


On  13  February  1988,  a  U.S.  Navy  (USN)  P-3C  Update  I  Orion 
aircraft  experienced  an  in-flight  failure  of  a  wing  leading- 
edge  section  during  a  high-speed ,  low-altitude  maneuver.  The 
aircraft  was  able  (with  some  difficulty)  to  return  to  its 
departure  point  and  make  a  safe  landing.  Some  three  years 
later,  on  27  April  1991,  an  Orion  operated  by  the  Royal 
Australian  Air  Force  (RAAF)  suffered  a  similar  but  more 
extensive  in-flight  failure  when  it  lost  three  wing  leading- 
edge  sections  during  a  maneuver  similar  to  the  U.S.  Navy 
mishap.  The  Australian  crew  was  unable  to  return  to  the 
runway  due  to  the  loss  of  lift  available  from  the  damaged 
wings.  The  aircraft  impacted  the  ocean  surface  short  of  its 
island  destination  in  a  nose  high-attitude  with  maximum  power 
set  on  all  four  engines.  One  crewmember  was  killed  when  a 
propeller  tore  loose  from  its  engine  and  entered  the  fuselage 
where  he  was  sitting. 

The  leading  edge  is  the  rounded  front  portion  of  the  wing 
which  is  physically  attached  to  the  front  wing  spar  and  is 
essential  in  the  production  of  lift  by  an  airfoil.  Loss  of  a 
.leading-edge  segment  results  in  a  dramatic  reduction  in  wing 
lift  capability  combined  with  a  corresponding  increase  in 


drag.  The  P-3  has  three  of  these  leading-edge  segments  on 
each  wing,  separated  by  the  engine  nacelles.  They  are  best 
referred  to  as  inboard,  center  and  outboard  sections.  Figure 
1  shows  the  P-3  aircraft.  The  inboard  sections  are  those 
between  the  fuselage  and  the 
inboard  engines  on  this 
four-engine  turboprop 
aircraft.  The  center 
sections  (blackened)  span 
the   distance   between   the 


inboard  and  outboard  engine 


Figure  1.  P-3  leading  edges. 


nacelles,  while  the  outboard 
sections  cover  the  remaining  distance  from  the  outboard 
engines  to  the  wing  tips.  The  length  (fore  and  aft)  of  these 
leading-edge  sections  is  15%  of  the  chord  (total  fore  and  aft 
distance)  of  the  wings. 

The  USN  P-3  lost  the  starboard  wing  center  leading-edge 
section  while  the  RAAF  aircraft  lost  both  of  its  center 
sections  and  the  starboard  inboard  section.  This  study 
focuses  on  the  center  leading  edge  sections. 

While  both  of  these  mishaps  were  investigated  by  the 
appropriate  authorities,  the  cause  of  the  failures  were 
undetermined,  though  widely  suspected  to  be  the  result  of 
aircraft  overstress  due  to  pilot  error.  In  both  cases,  the 
pilots  did  not  feel  that  they  had  exceeded  the  normal 
operating  envelope  (this  envelope  will  be  detailed  later)  for 


the  aircraft.  Even  if  the  aircraft  were  operated  outside  the 
authorized  envelope,  the  question  remains  as  to  why  the 
leading  edge  sections  failed  before  some  other  component. 

This  thesis  was  undertaken  for  the  purpose  of  studying  the 
aerodynamic  loading  and  structural  response  of  a  leading  edge 
section  due  to  operation  within  and  outside  the  normal  flight 
envelope.  Ultimately,  its  aim  was  to  determine  whether  there 
may  be  a  need  to  further  restrict  the  operating  envelope  or 
recommend  some  modification  to  the  existing  structure  in  order 
to  prevent  a  further  recurrence  of  the  in-flight  failure. 


II.   GENERAL  AERODYNAMICS 

A.   MISHAP  AIRCRAFT  SPECIFICATIONS 

The  configurations  and  operational  parameters  of  the 
mishap  aircraft  were  similar  in  that  both  were  at  high  gross 
weight  (RAAF  at  127,000  pounds,  USN  at  135,000  pounds)  and 
operating  at  high  airspeeds  (RAAF  approximately  380  knots,  USN 
approximately  350  knots).  Both  aircraft  executed  a  pullup 
maneuver  from  an  altitude  of  less  than  400  feet  at  the  time  of 
leading  edge  failure.  The  RAAF  maneuver  was  a  straight  pullup 
(wings  level),  while  the  USN  maneuver  was  a  starboard  rolling 
pullup  (meaning  that  the  aircraft  was  banked  to  the  right  as 
the  pullup  maneuver  was  executed) .  Interviews  with  some  USN 
crewmembers  indicate  that  the  failure  of  their  leading-edge 
section  may  have  begun  a  few  seconds  earlier  as  the  aircraft 
rolled  from  a  right  bank  to  wings  level  while  inbound  for  the 
above  stated  pullup.  This,  they  said,  was  reported  to  them  by 
ground  observers.  Table  1  presents  the  relevant  P-3C 
dimensional  data  used  during  this  analysis  while  Table  2 
delineates  the  aircraft  performance  parameters.  Figure  2 
shows  the  operating  flight  envelope  for  the  aircraft  in  a 
flaps-up,  gear-up  configuration. 


TABLE  1.   DIMENSIONAL  DATA  TRefs.  1,2,3 

Winer 

Area,  S  (ft2) 

1300 

Span,  b  (theoretical, 

ft)                      99 

MAC,  cbar  (ft) 

14.1 

Aspect  Ratio,  A 

7.5 

Taper  Ratio,  X 

0.4 

Dihedral,  (.25cw,  degrees)                    5.0 

Incidence,  Root  (degrees)                     3.0 

Tip 

0.5 

Airfoil  Section,  Root 

NACA0014-1.10  40/1.051  cli=.3,a= 

.8 

Tip 

NACA0014-1.10  40/1.051  cli=.4,a= 

.8 

Straight  Element 

0.15c 

Chord,  Root  (ft) 

18.9 

Tip 

7.6 

Aileron 

Area,  S^  (ft2) 

45.5 

Hinge  Line  (cw) 

0.725 

Deflection  Limit,  Up 

(degrees)               -23.3 

Down                       +16.2 

Horizontal  Tail 

Area,  St  (ft2) 

321.8 

Tail  Length  (.25cbarw 

to  .25cbart,  ft)          49.8 

TABLE  2.   PERFORMANCE  PARAMETERS  TRefs.  1,2,31 

Flight  Design  Gross  Weight  (pounds)  135000 

Load  Factor  (through  135,000  pounds,  G)  -1.0  to  +3.0 

Maximum  Operating  Speed,  Sea  Level  (knots)  405 

Center  of  Gravity  Limits  (135,000  lbs,  %MAC)  21.5  to  31.0 

Maximum  Shaft  Horsepower  (per  engine)  4600 

Maximum  Lift  Coefficient,  CL  x  (power  off)  1.30 

Lift  Curve  Slope,  CLa,  Tail  Off  (power  off)  4.84 

Tail  On  (power  off)  5.50 

Moment  Slope,  CmCL,  Tail  Off,  eg  @  0.25cbarw  0.19 
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Figure  2.  Operating  flight  envelope,  clean  configuration 
[Refc  2]. 


B.   PROBLEM  STATEMENT 

Accurate  analysis  of  the  problem  within  limited  time  and 
budget  constraints  available  for  completion  of  the  work  served 
to  restrict  the  options  for  solution  methods.  While 
instrumented  flight  and  laboratory  tests  would  probably  be 
more  precise,  these  choices  were  not  feasible.  It  was 
necessary  to  develop  an  analytic  method  which  could  provide 
credible  results  within  the  academic  environment. 

Determination  of  the  loads  applied  and  their  effect  upon 
the  leading  edge  structure  was  accomplished  by  a  three-step 
process.   The  first  step  was  to  conduct  a  static  aeroelastic 


span-load  analysis  of  the  wing  using  a  computer  program  to 
determine  the  section  lift  coefficients  and  structural  twist 
on  the  P-3  wing.  Next,  a  two-dimensional  panel  method  was 
employed  to  find  the  pressure  distribution  around  the  leading 
edge.  Finally,  the  forces  derived  from  the  pressure 
distribution  were  used  as  the  loads  applied  in  a  finite 
element  analysis  computer  application.  In  addition, 
consideration  was  given  to  the  aeroelastically-derived 
spanwise  wing  twists  which  were  introduced  as  torsion-like 
deflections  in  the  finite  element  analysis. 

C.   ASSUMPTIONS 

Certain  assumptions  were  made  in  conducting  the  analysis. 
While  static  aeroelasticity  effects  upon  the  wing  were 
accounted  for,  the  fuselage  and  tail  were  assumed  as  rigid 
structures.  In  addition,  the  effect  of  fuselage  interference 
on  wing  lift  distribution  was  neglected.  That  is,  the  wing 
was  taken  to  be  fully  effective  in  producing  lift  even  in  the 
central  region  where  it  is  influenced  by  the  fuselage.  The 
unswept,  straight-tapered  wing  with  an  aspect  ratio  of  7.5  was 
modeled  structurally  in  the  span-load  analysis  by  an  elastic 
axis.  Chordwise  bending  of  the  wing  was  neglected.  Inviscid 
solution  methods  were  employed  in  the  span-load  and  airfoil 
analysis  programs.  It  was  felt  that  these  assumptions  would 
model  the  flowfield  without  inducing  an  unacceptable  level  of 
error  in  the  overall  outcome,  and  that  this  analysis  was  a 


preferable  method  of  solution  to  performing  a  computational- 
fluid-dynamics  analysis  where  static  aeroelastic  influence 
could  not  be  included.  Static  aeroelastic  effects  upon  wing 
loadings  were  suspected  of  playing  a  large  role  in  the  problem 
because  of  the  location  of  the  wing's  elastic  axis  at  a 
constant  40  percent  of  chord,  according  to  available 
information.  This  is  well  aft  of  the  25  percent  of  chord 
location  at  which  a  majority  of  the  aerodynamic  loadings  are 
presumed  to  act,  creating  the  potential  for  a  significant 
coupling  influence  by  structural  deflection  during  the 
development  of  lift. 
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III.   WING  SPAN-LOAD  ANALYSIS 

A.   THEORY 

Because  the  P-3  wing  (as  well  as  nearly  all  others)  is 
flexible,  the  aerodynamic-load  distribution  varies  from  that 
which  would  be  seen  when  considering  the  wing  as  a  rigid 
structure.  Allowing  for  structural  deflection  of  the  wing 
increases  the  accuracy  of  the  result  in  the  attempt  to  model 
the  true  physics  of  the  problem.  The  computer  program  used 
for  this  analysis,  written  in  Microsoft  Quick  BASIC®,  was 
based  upon  the  work  of  Schmidt  [Ref.  4];  c.f.  Appendix  A.  The 
basic  concepts  employed  in  the  program  are  presented  below. 

Defining  a  set  of  linear,  simultaneous  equations  for  the 
span-load  solution  on  a  wing  at  a  specified  angle  of  attack 
starts  with  the  following: 

au(i/q)i  +  ai2(i/q)2  +.  .  .  +  ain(£/q)n  =  aL 
where  the  wing  is  cut  into  a  series  of  spanwise  stations  and 


•  a^  =  aerodynamic  influence  coefficient  of  ( j2/q )  at 

station  "j"  upon  induced  angle  at  control  station 
ii  -i  it 


•  2      =   running  span-load  (lb-in-1) 

•  q  =  freestream  dynamic  pressure  (lb-in"2) 

•  aL   =  geometric  angle  of  attack  at  control  station  "i" 


The  relationship  states  that  the  span-load-induced  downwash 
velocities  satisfy  flow  tangency  at  a  control  point. 
Expressed  as  a  matrix  equation  this  system  becomes: 

[A]{je/q}  =  {a} 
where, 

•  [A]  =  square  (n  x  n)  matrix  of  aerodynamic  influence 

coefficients  (length-1) 

•  {^/q}  =  column  (n  x  1)  matrix  of  span-load  values 

(length) 

•  {a}  =  column  (n  x  1)  matrix  of  angle-of -attack  input 

(radian) 

In  similar  fashion,  a  matrix  equation  relating  the 
structural  twist  at  a  wing  station  "i"  to  the  span-loading  may 
be  developed  from: 

sil£1   +  sL2£1   +  ...  +  sin£n   =   Aasi 
or,     q{sil(jg/q)1  +  si2(jg/q)2  +  ...  +  sin(£/q)n}  =  Aasi 
which  may  be  stated  to  apply  to  all  the  wing  stations  as 
before  to  yield  the  matrix  equation, 

q[S]{£/q}  =  {Aas} 
where, 


•  [S]  =  square  (n  x  n)  matrix  of  structural  influence 

coefficients  (length-force-1) 

•  {Aas}  =  column  (n  x  1)  matrix  of  effective  angle-of -attack 

changes  due  to  structural  twist  (radian) 


Next,  the  two  matrix  equations  may  be  combined  to  generate 
a  single  equation  for  an  elastic  wing  as  follows: 
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[A]{jg/q}E  =  {a}R  +{Aa}s 
[A]{je/q}E  =  {a}R  +  q[S]{£/q}E 
[  [A]  -  q[S]]  {Je/q}E  =  {a}R 
yielding  the  span-load  solution  for  the  elastic  wing, 

{je/q}E  =  [  [A]  -  q[S]]"1{a}R 
where  the  subscript  "E"  is  used  to  mean  elastic  and  "R"  to 
mean  rigid. 

The  concept  of  mathematical  symmetry  may  be  employed  to 
develop  equations  for  a  symmetric  and  anti-symmetric  load  case 
by  considering  the  wing  to  be  divided  into  left  and  right  hand 
wing  panels  at  the  spanwise  centerline.  This  process  allows 
a  superposition  of  linear  aerodynamic  solutions  such  as  the 
combination  of  a  symmetric  pullup  and  anti-symmetric  roll 
rate  to  yield  the  total  solution  for  a  rolling  pullup.  It 
also  reduces  computing  time  since  the  size  of  the  matrices  is 
half  the  original  required  for  the  total  wing.  In  employing 
this  technique  it  is  necessary  to  distinguish  between  the 
symmetric  and  anti-symmetric  forms  of  the  aerodynamic 
influence  coefficients. 

The  approach  employed  for  developing  the  aerodynamic 
influence  coefficients  involves  use  of  a  technique  known  as 
the  "Modified  Weissinger"  approach,  wherein  the  wing  is 
divided  into  a  series  of  spanwise  stations  with  swept  bound 
vortices  attached  at  the  local  quarter-chord  point,  giving 
rise  to  horseshoe  vortices  extending  downstream  to  infinity  in 
accordance  with  Helmholtz'  laws.   The  vortex  strengths  are 
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determined  in  accordance  with  the  Biot-Savart  law,  such  that 
the  induced  downwash  angle  at  all  local  three-quarter  chord 
points  (from  the  summation  of  all  vortices)  are  equivalent  to 
the  local  geometric  angle  of  attack.  The  three-quarter  chord 
points  are  termed  "control  points".  Enforcement  of  the  stated 
boundary  conditions  replicates  the  occurrence  of  tangential 
flow  over  the  surface  of  the  wing.  Application  of  the  Biot- 
Savart  and  Kutta-Jukowski  laws  to  the  geometric  relationships 
of  the  horseshoe  vortices  and  control  points  results  in 
development  of  the  symmetric  and  anti-symmetric  aerodynamic 
influence  coefficient  matrices, 

[As]  =  [  [K]m   +  [A]LH] 

'   [Aa3  =  [  [A]rh  "  [A] J 
where  the  subscripts  "RH"  and  "LH"  denote  the  right  and  left 

hand  wing  panels,  respectively.   These  influence  coefficient 

matrices  may  next  be  substituted  in  the  previously-developed 

equations  to  yield, 

[As]{£/q}    =  {ag} 

[Aa]{je/q}  =  {aa} 
Subsonic  compressibility  effects  were  included  in  the 
development  of  the  aerodynamic  influence  coefficients  using 
the  Prandtl-Glauert  planform  distortion  approach  [Ref.  5].  In 
this  method,  the  chordwise  dimension  of  the  planform  is 
increased  according  to  the  relationship, 
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x'=. 


\jl-M2 

The  "Modified  Weissenger"  approach  was  altered  to  a  panel 
form  for  this  analysis  in  that  the  wing  half  was  divided  into 
five   chordwise   and   ten   spanwise   stations,   requiring 
manipulation  of   (50  x 
50)    matrices.       A 
representative  sketch  of 
the  method  employed  is 
shown  in  Figure  3 .   The 
vector  U  in  the  figure 
represents    the    free 
stream  velocity  while  r 
is    the    circulation 


strength  of  the  vortex 
filament. 

Development  of  the 
structural  influence 
coefficient   matrix   in 


W\NG   DIVIDED 
INTO    P/VNELS 


HORSESHOE 
VORTEX   C1,J) 

PANEL   CON1ROL 
POINT 


Figure  3.  Wing  panel  model  [Ref  4] 


the  program  is  accomplished  through  determination  of  the 
moments  exerted  about  the  longitudinal  (X)  and  lateral  (Y) 
axes  of  the  wing  at  points  along  the  wing  elastic  axis  as 
shown  in  Figure  4.  These  points  correspond  to  the  mid-span 
locations  of  the  individual  panels,  along  the  elastic  axis. 
These  moments  are  found  by  solving  the  matrix  equations: 
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{My}=q(B/N)[t]{£/q} 
{Mx}=q(B/N)2[m]{jg/q} 
Where  B  is  the  span,  N 
is    the    number    of 
spanwise  local  stations, 
[t]   is   a   torsional 
matrix   and   [m]   is   a 
bending  moment  matrix. 
From  these  equations,  a 
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Figure  4.  Wing  moments  [Ref.  6] 


coordinate  axis  rotation  can  be  performed  as  follows: 

{T}  =  cosi\{My}  +  sin^{Mx} 
{M}  =  -  sini\,{My}  +  cosi\>{Mx} 
Next,  torsional  and  bending  stiffness  (GJ  and  EI)  values  along 
the  elastic  axis  are  employed  to  determine  the  angular 
deflections  resulting  from  the  torsion  and  bending  moments 
from  the  equations, 


y 

cosA, 


®t(Y)= 


I 


GJ(z) 


dz 


y 

cosA, 


et(m)= 


I 


M(z) 
EI(z) 


dz 


These  are  then  discretized  as, 


Aem(i)= 


B 


2/Vcos 


y>        M(J)    +   M(I) 


t  2NcosA^     ,.4^,    GJIJ)      GJ(I) 


j=i+1 


(J) 


(I) 
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or   expressed   in  matrix   form, 

{0t}= B    .    [u][  JL  ]{T\  {0J= B    .    [u]  [  JljM 

fc  ^iVcosA^  L  J  L  G J  J  "^  2NcosAe  L  J  L  EI  J 

The  angle-of -attack  column  vector  may  then  be  found  as, 

{&a3}   =  cos^{0t}  -  siniyem} 
which  is  also:    {Aas}  =  q[S]{l/q}. 

The  preceding  equations  may  then  be  tied  together  to  form  the 
equation  for  the  structural  influence  matrix, 

[S]  =  (B/2N)2[u][l/GJ][  cosi\>[t]  +  sin^(B/2N)  [m]J 
+  tani^(B/2N)2[u]  [l/EI]l  sini\>[t]  -  cosi\>(B/2N)  [m]J 

B.   USAGE  OF  THE  SPAN-LOAD  PROGRAM 

Application  of  the  wing  span-load  program  required  the 
determination  of  the  following  influences: 

•  Additional  loading  distribution  due  to  wing  angle  of 
attack 

•  Built-in  geometric  twist 

•  Airfoil  camber  distribution 

•  Dead-weight  induced  wing  twists  due  to  propulsion  system 
weight 

•  Aileron  float  angle 

•  Propeller  slipstream  effects 

•  Aileron  control  deflection 

•  Roll  helix  angle 
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The  last  two  influences  involved  ant i -symmetric  wing  span-load 
solutions.  The  total  wing  span-load  distribution,  which 
provided  a  measure  of  wing  lift  coefficient  (CL),  moment 
coefficient  (CmCG)  and  section  lift  coefficient  (C£) ,  was 
obtained  by  an  appropriate  linear  combination  of  the  above 
influences.  These  influences  were  incorporated  in  the  program 
as  discrete  angle  of  attack  adjustments  at  the  control  points. 
Finally,  a  specified  airplane  lift  coefficient  reguired  an 
estimate  of  the  tail  lift  contribution  before  the 
representative  wing  CL  could  be  estimated.  The  tail  lift 
contribution  was  based  upon  trimming  the  P-3  airplane  for  a 
specified  flight  condition  and  the  assumption  that  the  body 
and  horizontal  tail  behaved  approximately  as  rigid  structures. 
1.   Additional  Loading 

Static  aeroelastic  effects  upon  wing  span-loading  due 
to  geometric  angle  of  attack  were  incorporated  in  the  program 
as  a  selectable  input  from  the  operator  at  the  computer 
terminal.  Calculation  of  the  appropriate  angle  of  attack  for 
a  given  flight  condition  was  based  on  the  fundamental  equation 

CL  =  CL0  +  CLaa' 

with  appropriate  modification  for  tail  lift  contribution  as 
delineated  in  subsection  nine  of  this  Chapter. 

Early  in  the  analysis,  it  was  found  that  static 
aeroelastic  effects  had  a  dramatic  impact  on  additional 
loading.   Solutions  were  obtained  using  the  span-load  program 
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for  a  wide  range  of  dynamic  pressures.  The  q  equals  zero 
solution  corresponded  to  a  rigid-wing  case.  A  dynamic 
pressure  of  3.7  psi  corresponded  too  the  aircraft  operating  at 
a  Mach  number  (M)  of  0.6  at  sea  level.  The  variation  of  wing- 
alone  CLa  with  dynamic  pressure,  shown  in  Figure  5,  indicates 
a  41  percent  increase  due  to  static  aeroelastic  influences  at 
the  sea  level  flight  condition.  The  increase  in  lift-curve 
slope  is  associated  with  a  spanwise  variation  of  structural 
twist  as  shown  in  Figure  6.  At  a  Mach  number  of  0.6  for  sea- 
level  flight,  each  degree  of  geometric  angle-of -attack  input 
at  the  wing  root  results  in  1.95  degrees  of  a  at  the  tip,  with 
the  added  0.95  degrees  being  due  to  the  structural  twist 
component  in  the  a  direction.  The  static  aeroelastic 
influence  upon  the  wing  aerodynamic  center  was  determined  as 
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being  negligible,  a  result  which  may  be  attributed  to  the  wing 
0.25  chord  line  having  a  small  sweep  angle. 
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Figure  6.  Aeroelastic  variation  of  twist 
[Ref.  6] 


2 .   Built-in  Geometric  Twist 

The  effect  of  2.5  degrees  of  washout  was  included  in 
the  program  by  linearly  varying  the  local  (panel)  angle  of 
attack  moving  outward  from  zero  at  the  root  to  -2.5  degrees  at 
the  tip.  Figure  7  depicts  the  effect  of  washout  on  section 
lift  and  twist  distribution  for  the  rigid  and  elastic  P-3 
wing  cases.  The  magnitude  of  the  negative  section  C£  values 
is  seen,  in  Figure  7,  to  increase  by  21  percent  while  the  wing 
tip  experiences  a  1.5-degree  negative  twist  due  to  aeroelastic 
influences. 
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Figure  7.  Effect  of  built-in  geometric  twist  on  section 
lift  coefficient  and  twist  distribution. 

3.   Engine  Dead-weight  Twist 

A  separate  version  of  the  span-load  program  was 
developed  and  run  to  determine  the  wing  twist  induced  per  G 
due  to  the  effect  of  engine  dead-weight  moment.  Output 
information  was  generated  in  terms  of  discretized  angle-of- 
attack  adjustments  and  included  in  the  basic  span-load 
program.  An  engine  and  propeller  assembly  weight  of  3974 
pounds  was  assumed  to  act  at  the  (x,y)  coordinates  (coordinate 
system  as  shown  in  Figures  3  and  4)  of  (-51",  187")  and  (-45", 
357")  for  the  inboard  and  outboard  engines,  respectively. 
Figure  8  shows  the  effect  of  engine  dead-weight  twist  on  the 
elastic  twist  distribution  at  three  G's  and  405  knots.   The 
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wing  is  twisted  to  a  maximum  of  -4.5  degrees  at  the  tip  under 
this  load  condition. 
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Figure  8.  Effect  of  engine  dead-weight  on  structural 
twist . 

4 .   Camber 

Although  listed  in  Table  1  as  a  NACA  0014  which  tapers 
to  a  NACA  0012,  the  airfoil  used  for  the  P-3  wing  is  not  a 
symmetric  section,  as  indicated  by  available  data  [Ref .  8,  p. 
2-470].  It  is,  in  fact,  a  hybrid  with  considerable  camber. 
The  wing  camber  has  an  effect  on  aeroelastic  behavior  and  was 
included  in  the  analysis  by  determination  of  camber-line  slope 
at  the  five  chordwise  control  points  which  make  up  the  wing 
model  and  entering  the  negative  of  this  value  as  an  adjustment 
to  the  station  angle  of  attack. 
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5.  Aileron  Float 

The  span-load  analysis  included  the  effect  of  aileron 
float  as  given  in  a  graph  of  aileron  angle  versus  airspeed  in 
a  Lockheed  report  [Ref.  9,  p.  123].  An  equation  curve  fit  for 
the  50-pound  control-wheel-force  curve  was  developed  and 
incorporated  in  the  program.  Deflected  aileron  surface  area 
was  matched  by  a  selection  of  panels  on  the  wing  model  for 
varying  degrees  of  deflection  in  order  to  closely  emulate  the 
aircraft  control-surface  deflection  and  aeroelastic  effect. 

6.  Aileron  Deflection  Angle 

Using  the  same  process  as  that  used  for  the  aileron 
float  angle,  the  effect  of  aileron  deflection  for 
consideration  of  anti-symmetric  solutions  to  emulate  roll 
maneuvering  was  also  included.  A  curve  fit  to  the  data  given 
for  available  50-pound  control-wheel-force  deflection  in  the 
same  report  [Ref.  9,  p.  123]  was  used  to  generate  the 
deflection  angles. 

7.  Roll  Helix  Angle 

Data  [Ref.  9,  p.  118]  for  available  tip  helix-angle 
(pb/2V)  variation  with  airspeed  were  curve  fitted  and 
incorporated  as  another  control  point  angle-of-attack 
variation  in  the  program.  A  roll  helix  angle  of  2.63  degrees 
provided  roll  moment  equilibrium  with  the  available  9.3 
degrees  of  aileron  deflection  at  275  knots  equivalent 
airspeed,  at  a  control  wheel  force  of  50  pounds.   At  405 
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knots,  the  roll  helix  angle  was  1.03  degrees  for  roll 
equilibrium  with  an  aileron  deflection  of  8.2  degrees.  The 
merging  of  the  anti-symmetric  input  due  to  aileron  deflection 
with  the  symmetric  pullup  solution  allowed  a  comprehensive 
analysis  of  the  static  aeroelastic  effect  of  a  rolling  pullout 
maneuver  as  described  in  part  A  of  this  Chapter.  An  example 
of  the  individual  and  combined  effects  of  these  components  on 
wing  section  lift  coefficient  (Cjg)  distribution  is  shown  in 
Figure  9,  where  the  outer  corner  of  the  operating  envelope  for 
roll  maneuvering  (405  knots,  2.4  G's)  was  examined. 
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Figure  9.  Rolling  pullup  contributions  to  section  lift 
distribution . 

8.   Propeller  Slipstream  Effect 

The  increase  in  dynamic  pressure  over  the  wing  due  to 
the  propellers  was  calculated  using  momentum  theory  as  stated 
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in  Glauert  [Ref.  10,  p.  200]  according  to  the  following 
formula: 

T  =  Ap(V  +  v-^v-l 
where, 

•  T  =  propeller  thrust  (pounds) 

•  A  =  propeller  disk  area  (ft2) 

•  V  =  freestream  velocity  (ft/sec) 

•  vx  =  velocity  increase  behind  the  propeller,  determined 

from  known  thrust. 

The  velocity  increase,  vlf  was  converted  to  a  dynamic  pressure 
boost  and  incorporated  in  the  span-load  program  in  the  area 
behind  the  propellers. 

9 .   Tail  Contribution 

Horizontal  tail  and  fuselage  moment  effects  on  wing 
span-load  distribution  resulted  in  an  adjustment  to  the  wing 
angle-of -attack  value  used  as  input  to  the  program.  This 
adjustment  was  calculated  according  to  the  relationship, 


CL  ~  CLadd  +  CLt  +  CL0 


where, 


•  CL  =  lift  coefficient  required  for  flight  condition 

•  CLadd  =  lift  coefficient  due  to  additional  loading  (due  to 

angle  of  attack) 

•  CLt  =  tail  contribution  to  lift  coefficient 

•  CL0  =  tail-off  lift  coefficient  at  zero  angle  of  attack 

(due  to  camber,  wing  twist  and  dead-weight  twist) 
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and, 

CLt  -  -ACmt(cbar/i>t) 
where, 


•  ACmt  =  tail-moment  coefficient 

•  cbar  =  mean  aerodynamic  chord  of  the  wing 

•  J0t  =  distance  from  .25cbarw  to  . 25cbart 


Additionally,  at  airframe  trim  (CmCG  =  0), 


Cmt    Cm0  +  CmCLCLadd 


where, 


•  Cm0  =  moment  coefficient  at  zero  angle  of  attack 

•  CmCL  =  tail-off  airplane  dCm/dCL,  e.g.  =  0.25cbar 


The  above  was  assembled  and  solved  for  CLadd  to  give, 


Ladd' 


tCL  CL0]     — 7 Cm0 


1  + 


cbar 


mCL 


The  value  of  tail-off  CmCL  for  the  aircraft  was  found  from 
wind  tunnel  data  to  be  approximately  0.19  [Ref.  3,  p.  20]. 
The  value  of  CLadd  for  any  given  flight  condition  was  then  used 
to  solve  for  the  angle  of  attack  in  the  span-load  program 
according  to  the  formula: 

a  =  CLadd/CLa 
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The  value  of  CLa  was  available  as  an  output,  for  the  elastic 
wing,  from  the  program  for  any  given  airspeed  by  entering  one 
radian  for  the  angle  of  attack  input. 

A  correction  of  +0.03  to  Cm0  due  to  fuselage  effect 
was  made  after  finding  that  the  value  Cm0  in  the  wind  tunnel 
data  and  that  from  the  program  differed  by  this  amount.  This 
correction  was  verified  by  calculations  made  in  accordance 
with  Etkin  [Ref.  7 ,  p.  334]  concerning  the  effect  of  body  and 
engine  nacelles  on  neutral  point  location. 

C.   EFFECTS  NOT  INCLUDED 

Two  other  factors  were  considered  as  possible  contributors 
to  the  static  aeroelastic  problem,  but  were  found  to  have  no 
significant  impact  and  therefore  not  included  in  the  span-load 
program.  These  factors  were  propeller  gyroscopic  effects  and 
the  effect  of  moments  due  to  wing  fuel. 
1.   Propeller  Gyroscopic  Moment 

This  phenomenon  was  investigated  using  the  following 
formulation  from  [Ref.  11]: 

M  =  2Ipojfi 
Ip  =  mk2 
where, 

•  M  =  moment  due  to  gyroscopic  precession  (lb-ft) 

•  Ip  =  polar  moment  of  inertia  of  each  blade  (lb-ft  sec2) 

•  (0  =  rotation  rate  of  the  propeller  (rad-sec-1) 
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•  Q  =  pitch  rate  of  the  aircraft  (rad-sec-1) 

•  m  =  mass  of  propeller  blade  (slugs) 

•  k  =  radius  of  gyration  (ft) 

The  analysis  showed  that,  allowing  for  a  five-degree-per- 
second  aircraft  pitch  rate,  the  moment  developed  by  each 
propeller  was  1225  lb-ft  in  a  counterclockwise  direction  as 
viewed  from  above  the  propeller.  This  moment  was  considered 
insignificant  since  it  is  applied  in  the  lateral  plane  and 
does  not  influence  the  static  aeroelastic  span-load  solution. 
2  c   Wing  Fuel  Moment 

Analysis  of  the  fuel  tank  geometry  and  location 
revealed  that  the  center  of  mass  of  the  fuel  (with  full  wing 
tanks)  lies  nearly  coincident  with  the  elastic  axis.  Any 
torsional  moment  derived  from  this  source  would  be  negligible. 

D.   APPLICATION 

Linear  summation  of  the  contributing  factors  addressed 
earlier  in  this  Chapter  resulted  in  a  prediction  of  the  total 
span-load  distribution  for  the  P-3  wing  under  any  given  flight 
condition.  The  primary  flight  conditions  of  concern  in  this 
analysis  were  those  encountered  during  the  aircraft  mishaps  as 
discussed  in  the  introduction.  The  basic  premise  applied  was 
to  examine  limit  loads  at  the  edge  of  the  operating  envelope 
and  then  expand  the  analysis  to  regions  outside  the  envelope, 
at  the  ultimate  load  condition.  In  addition  it  was  decided  to 
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first  examine  the  loads  encountered  in  a  symmetric  pullup  of 
three  G's  at  275  knots.  This  speed  was  chosen  because  of  an 
observed  higher  load  on  the  leading  edge  than  determined  at 
the  upper  right  corner  of  the  flight  envelope,  at  405  knots. 
(A  further  discussion  of  this  loading  phenomenon  will  be 
presented  in  Chapter  IV. )  Following  that,  the  target  flight 
condition  was  extended  to  4.5  G's  at  325  knots,  which 
represented  an  approximate  extension  of  the  envelope  using  a 
value  for  CLmax  of  1.3.  Next,  the  effect  of  anti-symmetric 
loading  in  the  form  of  a  starboard  and  then  a  port  rolling 
pullup  were  considered  in  order  to  assess  the  contributions  of 
aileron  deflection  and  roll  helix  angle.  Rolling  pullup 
analyses  were  done  at  2.4  G's  and  275  knots  to  remain  inside 
the  envelope  and  maintain  some  congruity  with  the  symmetric 
loading  case.  Presented  in  this  section  are  plots  of  the 
impact  of  some  of  the  various  flight  maneuvers  on  section  lift 
coefficient  and  twist  distribution,  using  the  fully-developed, 
tail-on  solution  with  all  contributing  factors  included. 
Figures  10  and  11  show  the  275-knot,  3-G  symmetric  condition, 
and  compare  section  lift  coefficient  and  structural  twist  for 
the  rigid-wing  and  elastic-wing  cases.  In  Figures  12  and  13, 
the  same  distributions  are  depicted  for  the  275-knot,  2.4-G 
rolling  pullup  load  condition. 
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Figure  10.  Lift  coefficient  distribution  at  275  knots, 
3-G,  symmetric  pullup. 
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Figure  11.  Twist  distribution  at  275  knots,  3-G, 
symmetric  pullup. 
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Figure  12.   Lift  coefficient  distribution  at  275  knots, 
2.4-G,  rolling  pullup. 
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Figure  13.  Twist  distribution  at  275  knots,  2.4-G, 
rolling  pullup. 
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IV.   TWO-DIMENSIONAL  PRESSURE  ANALYSIS 

A.   THEORY 

Like  the  static  aeroelastic  span-load  analysis,  the  two- 
dimensional  pressure  analysis  employed  in  this  thesis  involves 
the  application  of  linear  superposition.  This  is  a 
consequence  of  the  pressure  analysis  being  based  upon  a 
solution  to  the  LaPlace  equation,  a  linear,  homogeneous 
second-order  partial  differential  equation.  Linearity  allows 
the  problem  to  be  subdivided  into  three  separate  elements 
(which  will  be  described  later  in  this  Chapter)  and  added 
together.  The  actual  application  was  based  upon  a  panel 
method  similar  to  the  technique  used  in  the  preceding  Chapter 
except  that  now,  instead  of  dividing  a  wing  planform  into 
chordwise  panels,  a  two  dimensional  (x,y)  airfoil  is  divided 
into  panels  along  the  perimeter  of  its  surface.  Camber  is 
defined  in  the  airfoil  shape,  instead  of  being  added  on  as  a 
discretized  variation  of  angle  of  attack,  as  before.  The 
analysis  now  concerns  a  vertical  plane  or  cross-section. 
Steady  (no  variation  in  the  flow  field  with  time),  inviscid, 
incompressible  flow  is  assumed  to  exist.  The  panel  method  and 
formulation  employed  is  documented  in  a  Naval  Postgraduate 
School  thesis  by  Teng  [Ref.  12]. 
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1 .  Coordinate  Axis  and  Panel  Numbering  System 

The  airfoil  is  considered  to  be  fixed  in  an  (x,y) 
coordinate  system  with  its  origin  at  the  point  of  intersection 
of  the  chord  line  with  the  leading  edge.  The  positive  x  axis 
points  aft  to  the  trailing  edge  while  positive  y  is  up.  The 
panels  which  make  up  the  airfoil  surface  are  of  varying 
lengths,  depending  mostly  on  the  radius  of  curvature,  and  are 
numbered  from  1  through  n  starting  at  the  lower  surface  of  the 
trailing  edge  and  proceeding  clockwise  to  the  upper  surface  at 
the  trailing  edge.  Delineating  the  end  points  of  these  panels 
are  nodes  which  begin  with  the  number  1  at  the  trailing  edge 
and  proceed  along  the  same  numbering  path  as  the  panels.  The 
trailing  edge  point  is  counted  twice,  giving  n+1  nodes  in  all. 

2 .  Flow  Formulation 

Consider  some  panel  j  on  the  surface  of  the  airfoil. 
On  this  surface  there  exists  a  pair  of  singularity 
distributions,  known  as  a  source  distribution  q^  and  a 
vorticity  distribution  y.  The  strength  of  the  source 

distribution  varies  from  panel  to  panel  while  the  vortex 
strength  is  the  same  for  all  panels.  These  singularity 
distributions  satisfy  LaPlace's  equation  and  the  far  field 
boundary  condition. 

Applying  superposition,  the  overall  flow  field  is 
considered  to  be  made  up  of  three  individual  flows  and  is 
represented  by  the  equation, 
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where  0,,,  is  the  potential  of  the  freestream  flow, 

0oo  =  voo(x  cosa  +  y  sina) 
<t>3    is  the  velocity  potential  of  the  source  distribution  of 
strength  q(s)  per  unit  length  (s)  and  is  calculated  by, 


4>s-j^lln(r)ds 


where  (r)  is  the  radial  distance  from  some  point  at  which  a 
source  and  vortex  flow  exist,  to  the  midpoint  of  the  panel  in 
consideration.  In  addition,  <pv  is  the  velocity  potential  of 
a  vorticity  distribution  of  strength  A(s)  per  unit  length  and 
is  given  by, 


*v~m£**s 


where  0  is  the  angle  formed  by  a  line  drawn  along  the  radial 
distance  (r)  and  the  panel  in  question. 

Each  of  the  preceding  equations  is  integrated  along 
the  straight  line  which  makes  up  each  of  the  panels,  where  qj 
and  X are  constant.  The  individual  effects  are  then  summed  to 
give  the  total  effect  of  the  sources  and  vortices  from  all 
panels,  as  given  in  the  equation, 

n 
$=V0D(xcosa 


+ysina)  +  £       [  ||ln(r)  --^0]ds 

J=l  panel  ( j) 
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The  calculation  of  $  requires  solution  of  the  (n+1) 
unknowns,  qj  (j  =  1,2, ...,n)  and  y«  This  is  accomplished 
numerically  in  the  computer  program.  Once  $  is  known,  the 
velocity  can  be  found  by  taking  the  gradient  (V)  of  $.  The 
total  velocity  vector  is  found  as, 

Vtotal  "  V*  "  V<*>»  +  V(0S  +  0V) 
Next,  the  coefficient  of  pressure  is  found  from  the 
Bernoulli  equation  in  the  incompressible  form, 

/-.  =1  _  /  " total  N  2 

3 .   Boundary  Conditions 

Both  the  condition  of  flow  tangency  at  the  surface  and 
the  Kutta  trailing  edge  condition  must  be  satisfied  as 
boundary  conditions.  As  in  the  span-load  analysis,  control 
points  are  designated  at  which  flow  tangency  is  satisfied, 
except  here  the  control  point  is  taken  as  the  mid  point  of  the 
panel.  It  is  stipulated  that  each  control  point  will  have 
tangential  velocity,  (Vt)i,  but  that  all  normal  velocities, 
(Vn)i,  will  be  exactly  zero. 

The  Kutta  condition  requires  that  the  pressures  on  the 
upper  and  lower  surface  at  the  trailing  edge  be  equal.  Using 
Bernoulli's  equation  for  steady  potential  flow,  this  state  of 
pressure  equilibrium  is  found  to  exist  when  the  tangential 
velocities  in  the  downstream  direction  are  equal  at  the  upper 
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and  lower  trailing-edge  panels.  In  equation  form  this  is 
written, 

(Vt)1   =  ~(Vt)n 
The  task  then  becomes  one  of  using  the  boundary  conditions  to 
solve  for  the  (n+1)  unknowns. 
4.   Influence  Coefficients 

The  concept  of  influence  coefficients  is  again 
employed  in  this  portion  of  the  analysis,  as  it  was  in  the 
span-load  analysis.  Here,  the  influence  coefficients  take 
the  form  of  induced  normal  and  tangential  velocities  at  the 
control  point  of  a  given  panel.  These  velocities  are  induced 
by  the  source  and  vorticity  distributions  of  the  other  panels, 
and  it  is  from  this  influence  that  they  receive  their 
designations: 

•  A":  j  =  normal  velocity  induced  at  the  ith  panel  control 
point  by  the  source  distribution  on  the  2       panel. 

•  At;j  j  =  tangential  velocity  induced  at  the  ith  panel  control 
point  by  the  source  distribution  on  the  jth  panel. 

•  B1^ j  =  normal  velocity  induced  at  the  ith  panel  control 
point  by  the  vorticity  distribution  on  the  jth  panel. 

•  Bti  j  =  tangential  velocity  induced  at  the  ith  panel  control 
point  by  the  vorticity  distribution  on  the  jth  panel. 

These  influence  coefficients  are  calculated  through 
application  of  the  geometric  relationships  which  exist  between 
the  panels  in  conjunction  with  the  formulas  for  the  source  and 
vorticity  velocity  potentials  as  given  above. 


34 


5.   Numerical  Solution  Method 

Using  the  influence  coefficients,  the  boundary 
conditions  can  now  be  employed  to  write  (n+1)  equations  which 
may  be  solved  in  matrix  form  for  the  (n+1)  unknowns.  The  set 
of  n  equations  comes  from  enforcement  of  the  flow  tangency 
boundary  condition  in  the  form, 

n  n 

£  [A'i^lnE  Bn±j+Vws±n(a-Q±)=0 
j=l  j=l 

Next  may  be  written  the  enforcement  of  the  Kutta  boundary 
condition  as, 

n  n  n  n 

"E  [A'ij^-yE  Btlj-7<ncos(a-01)=X;  [AtnjgJ.]+YX)B%'+^cos(a"ei) 
j-1  j=l  j=l  j=l 

The  negative  signs  on  the  left  side  of  the  equation  are  due  to 
the  defined  orientation  of  the  tangential  velocities  as 
positive  in  the  downstream  direction. 

These  equations  may  then  be  expressed  in  matrix  form 
with  the  (n+1)  unknowns  (i.e.,  qj  ( j=l, 2, . . . ,n)  and  y) 
arranged  as  a  column  (nxl)  matrix  multiplied  by  the  An+1  n+1 
(n+1  x  n+1)  influence  coefficient  matrix  and  set  equal  to  a 
Bn+1  column  matrix.  From  this  point,  a  Gaussian  Elimination 
numerical  technique  may  be  employed  to  solve  for  the  (n+1) 
unknowns  [Ref.  13]. 
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6.   Velocity  and  Pressure  Distribution 

Once  the  q^  and  y  are  found,  the  tangential  velocities 
at  the  control  points  can  be  solved  for  according  to  the 
equation: 

n  n 

where  i  =  l,2,...,n.    From  this,  the  individual  pressure 
coefficients  may  be  solved  using  the  equation, 

(Cp)i  =  1  -  (V1^2,   i=l,2,...,n 
At  this  point,  the  forces  at  work  on  the  airfoil  may  be  found 
by  first  integrating  forces  in  the  airfoil  coordinate  system 
as  follows, 

n 
cy=_E  (Cp)±(xi+1-x±) 
i=l 

n 

cx=E  (Cp)i(yi+i-yi) 

i=l 

Performing  a  coordinate  axis  rotation  to  re-align  with  that  of 
the  freestream  yields  the  lift  coefficient, 

C£  =  C  cose?  -  Cx  sina 
For  this  application,  the  values  of  Cx  an  C  are  found  in 
discretized  form  at  each  panel  as  (Cx)i  and  (C  )ir  i  = 
l,2,...,n,  and  then  multiplied  by  the  chord  length  and 
freestream  dynamic  pressure  to  arrive  at  the  normal  and 
chordwise  force  exerted  at  each  panel.    From  this,   the 
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leading-edge  panels  are  selected  and  their  forces  collected 
for  application  as  point  loads  in  the  finite  element  analysis 
portion. 

B.   APPLICATION 

1.  P— 3  Airfoil  Section 

As  mentioned  in  Chapter  III,  the  airfoil  shape  used  in 
this  application  is  as  delineated  by  [Ref.  8].  Surface 
coordinate  locations  were  solved  using  the  tables  and 
equations  provided  therein  relative  to  a  "wing  reference 
plane"  which  appeared  to  correspond  to  a  water  line  (i.e.  the 
angle  of  incidence  at  the  root  was  included  in  the 
definition).  These  coordinates  were  then  rotated  to  an  (x,y) 
coordinate  system  aligned  with  the  chord  of  the  airfoil  as 
required  in  this  panel  method.  The  result  was  an  airfoil 
consisting  of  some  48  panels,  which  was  increased  to  52  panels 
(53  node  points)  due  to  observed  roughness  of  the  leading  edge 
shape  when  plotted.  This  smoothing  of  the  leading  edge  shape 
was  achieved  by  applying  the  specified  leading  edge  radius  to 
create  intermediate  node  points.  The  basic  shape  of  the 
airfoil,  though  not  precisely  to  scale,  is  depicted  in  Figure 
14.   The  locations  of  the  node  points  are  also  shown. 

2 .  Program  Inputs  and  Outputs 

The  desired  output  from  the  two-dimensional  panel 
method  program  was  a  collection  of  forces  and  boundary 
constraints  to  be  applied  at  node  locations  in  the  finite 
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Figure  14.  P-3  airfoil  showing  panel  nodes.  (Not  to 
scale) 

element  analysis  phase.  In  order  to  achieve  this,  the  program 
was  altered  considerably  from  the  initial  state  as  outlined 
earlier  in  this  Chapter.  The  final  FORTRAN  code  is  available 
in  Appendix  B.  Since  it  was  necessary  to  develop  loads  to  be 
applied  to  a  three-dimensional  model,  the  program  was  set  up 
to  iteratively  compute  the  force  distribution  at  a  series  of 
two-dimensional  airfoil  sections  which  ranged  in  location  from 
the  inboard  end  of  the  finite  element  model  at  wing  station 
(WS)  256  to  the  outboard  end  at  WS  320.  (These  locations 
correspond  to  the  outboard  64  inches  of  the  wing  center 
section  leading  edge,  located  between  the  nacelles.  )  This 
section  force  distribution  was  scaled  up  from  a  chord- 
normalized  airfoil  shape  to  a  full-sized  airfoil  as  described 
above,  then  multiplied  by  a  scaling  factor  equal  to  the 
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distance  between  the  section  locations.  Section  locations 
consisted  of  24  wing-station  positions  along  the  leading-edge 
model,  including  the  nine  ribs  and  15  intermediate  points  as 
determined  by  finite  element  model  node  locations. 

The  first  of  the  inputs  to  the  program  consisted  of 
the  airfoil  geometry  definition  as  described  above.  Next,  a 
tabular  file  of  wing  station  locations  was  read  in  together 
with. the  section  lift  coefficient  for  that  spanwise  location 
as  found  by  the  span-load  program.  These  section  lift 
coefficients,  initially  found  for  the  .35,  .45  and  .55t| 
locations  (where  T|  =  2y/b)  were  curve-fitted  using  the  Cricket 
Graph  plotting  software  in  order  to  achieve  a  high  degree  of 
accuracy  in  determining  the  individual  Cjg's  at  stations  which 
were  no  more  than  three  inches  apart.  Also  included  in  this 
input  file  was  a  list  of  spanwise  multiplication  factors  for 
scaling  up  the  load  as  described  above.  Additional  input 
files  consisted  of  the  finite  element  node  numbers  which  were 
matched  with  their  respective  x  and  y  direction  loads  in  the 
program.  Read  in  from  the  terminal  were  the  airspeed  under 
consideration  and  the  twist  angles  of  the  wing  box  as 
determined  from  the  span-load  program. 

Output  consisted  primarily  of  the  load  file  which  included 
not  only  the  loads  at  the  finite  element  node  points,  but  also 
the  constraints  and  twist  displacements  for  the  upper  flange 
and  lower  hinge  node  points  of  the  leading-edge  finite  element 
model.   The  twist  displacements  were  calculated  within  the 
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panel  program  using  formulas  based  on  geometric  considerations 
of  the  height  of  the  wing  spar  at  the  inboard  and  outboard 
ends  of  the  leading-edge  segment,  and  the  net  twist 
displacement  of  the  front  spar  from  the  inboard  to  the 
outboard  end  of  the  model.  That  is,  the  finite  element  model 
was  assumed  to  have  undergone  a  rigid-body  rotation  to  the 
degree  of  twist  which  was  found  to  exist  at  the  inboard  end 
(WS  256).  Therefore,  twist  displacements  at  the  inboard  end 
were  set  to  zero,  followed  by  application  of  the  subsequent 
net  twist  distribution  that  occurred  at  the  other  wing 
stations  while  proceeding  outboard  to  WS  320.  This  net  twist 
distribution  was  equal  to  the  difference  between  the  outboard 
and  inboard  twist  amounts,  applied  linearly  over  the  24 
stations.  This  twist  amounted  to  approximately  0.3  degrees  in 
the  275-knot,  3-G  symmetric  pullup.  The  displacements 
generated  for  application  at  the  node  points  were  in  the 
longitudinal  (x)  and  vertical  (z)  directions  on  the  three- 
dimensional  model.  Twist  rotation  of  the  front  spar  was  taken 
to  be  about  its  vertical  mid  point,  and  the  structural  wing 
box  was  assumed  to  have  no  chordwise  distortion  as  it  rotated 
about  the  .40c  elastic  axis  location.  Other  output  took  the 
form  of  files  to  examine  C  distributions  and  to  tabulate 
loads  by  wing  station  and  two-dimensional  node  point  for 
verification  of  the  load  file.  In  addition,  output  was 
generated  which  approximated  the  total  normal  and  chordwise 
loads  applied  to  the  entire  wing  leading-edge  center  section 
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located  between  the  engine  nacelles.    This  estimation  was 
accomplished  by  computing  the  load  on  the  leading-edge  portion 
of  the  airfoil  at  WS  274  and  multiplying  by  the  length  of  the 
leading-edge  segment  (92  inches). 
3 .   Program  Operation 

After  reading  in  the  coordinate  information  for  the 
airfoil  along  with  the  wing  station,  C£  and  load 
multiplication  factor,  the  FORTRAN  program  calculated  the 
chord  length  at  the  particular  station  based  upon  the  0.40 
taper  ratio.  It  also  determined  the  thickness  fraction  of  the 
airfoil  section  at  that  station  by  assuming  a  linear  taper 
from  14%  maximum  thickness  at  the  root  to  12%  at  the  wing  tip. 
This  thickness  factor  was  then  used  to  recalculate  the  y 
coordinate  position  of  each  node  point  to  redefine  the  shape 
of  the  airfoil.  An  initial  angle  of  attack  of  one  degree  was 
set  and  the  process  described  in  the  theory  section  of  this 
Chapter  took  place,  wherein  the  Cx,  C  ,  and  C£  were  calculated 
for  that  angle  of  attack.  This  C£  value  was  then  compared  to 
that  required  (as  input  with  the  wing  station),  and  an 
iterative  cycle  commenced  in  which  the  angle  of  attack  was 
varied  up  or  down  by  an  amount  based  on  the  product  of  the  C2 
deviation  multiplied  by  a  preset  angular  value.  An  accuracy 
test  of  .0001  was  applied  to  reach  an  acceptable  value  for  C£, 
at  which  time  the  process  started  over  with  the  next  wing 
station.   Forces  in  the  x  and  y  direction  were  matched  with 
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the  appropriate  finite  element  node  points  and  written  to  the 
load  file  for  each  iterative  cycle.  After  all  loads  were 
calculated  and  stored,  the  twist  displacements  and  zero 
boundary  constraints  were  calculated  and  appended  to  the  load 
file. 

4.  Verification  of  the  Program 

The  accuracy  of  the  two-dimensional  panel  method  was 
verified  by  comparison  with  published  empirical  data  for 
tangential  velocity  and/or  pressure  coefficient  distributions 
for  the  NACA  0012  and  Eppler  E64  airfoils  before  its  use  in 
this  application.  Results  were  nearly  identical  to  the 
published  data  with  only  a  small  deviation  seen  near  the 
trailing  edge  of  the  program  tangential  velocity  distribution 
for  the  Eppler  airfoil.  No  difference  from  the  NACA  0012  Cp 
data  could  be  identified. 

5.  Flight  Regime  Selection 

It  was  found  in  the  course  of  running  the  program  at 
various  airspeeds  and  angles  of  attack  that  the  highest  loads 
on  the  leading-edge  segment  were  generated  at  slower  airspeeds 
and  higher  angles  of  attack  as  a  constant  G  load  was 
maintained  on  the  aircraft.  This  result  seemed  contrary  to 
conventional  opinion  that  the  highest  loads  would  most  likely 
occur  at  or  near  the  high  speed  end  of  the  operating  envelope. 
A  brief  study  was  undertaken  to  determine  the  cause  of  this 
phenomenon  and  a  hypothesis  is  given  here. 
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a.      Method  Employed 

In  an  effort  to  study  the  effect  of  dynamic 
pressure  (airspeed)  and  angle  of  attack  on  the  leading  edge  of 
an  airfoil,  the  variation  of  net  Cp  distribution  on  the 
leading  edge  at  various  angles  of  attack  was  first  examined. 
These  Cp  values  were  the  numerical  sum  of  the  difference 
between  the  upper  and  lower  Cp's,  integrated  over  the 
chordwise  distances  occupied  by  their  respective  panels. 
These  data  obtained  were  then  curve  fitted  with  a  third  order 
polynomial  and  used  in  a  spreadsheet  to  calculate  the  Cp 
distribution  at  varying  angles  of  attack.  These  angles  of 
attack  were  generated  by  varying  the  airspeed  from  275  to  425 
knots,  calculating  dynamic  pressure  (q)  for  a  3-G  wing  loading 
from  the  Bernoulli  equation,  and  then  converting  these  q's  to 
angles  of  attack  required  using  the  basic  lift  formula, 
altered  by  the  equation, 

CL  =  CL0  +  CLaa 

to  give, 

a  =  3W/CLa<3S    and    «  =  3W/CLa<3S  "  CL0/CLa 

for  the  symmetrical  NACA  0012  and  cambered  P-3  airfoil, 
respectively.  These  angles  of  attack  were  then  used  as  inputs 
to  the  polynomial  curve  fits,  from  which  a  corresponding  set 
of  C  values  were  calculated.  Next  the  product  of  Cp  and  q 
were  found,  to  give  the  pressure  acting  on  the  leading  edge, 
corresponding  to  a  matched  set  of  q  and  angle  of  attack.  This 
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information  was  then  plotted  as  seen  below  in  Figures  15  and 
16. 

jb.  Symmetric  Airfoil 

Figure  15  shows  the  distribution  of  the  various 
parameters  on  a  relative  scale  as  airspeed  increases  for  the 
case  of  a  symmetric  airfoil.  Note  that  the  pressure  acting  on 
the  leading  edge  is  nearly  constant,  showing  only  a  slight 
increase  with  increasing  airspeed. 


Pressure  Change  of  0012  Airfoil  LE  with  AOA  and  Q 
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Figure  15.  Pressure  analysis  of  0012  leading  edge. 

c.  Cambered   P-3  Airfoil 

Next  the  same  information  is  plotted  for  the  P-3 
airfoil  in  Figure  16.  Note  that  the  pressure  on  the  leading 
edge  undergoes  a  marked  decline  as  airspeed  increases  (AOA 
decreases) . 
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Pressure  Change  of  P-3  Airfoil  LE  with  AOA  and  Q 
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Figure   16.    Pressure   analysis  of   P-3    leading  edge. 

d.      Analysis 

Because  of  the  apparent  drop  in  pressure  on  the 
cambered  airfoil,  it  was  concluded  that  the  effect  of  camber 
was  to  shift  the  loading  of  the  wing  in  a  way  that  caused 
angle  of  attack  to  become  the  dominant  influence  rather  than 
dynamic  pressure.  It  was  also  observed  that  if  the  lift 
eguation  applied  in  the  program  was  that  of  a  cambered  airfoil 
(i.e.  Cjg0  had  some  positive  value  as  the  lift  curve  was 
displaced  upward),  this  loading  phenomenon  was  present.  If 
C^q  were  zero,  the  load  on  the  leading  edge  was  approximately 
the  same  at  various  angles  of  attack  and  airspeeds.  Because 
of  this  observation,  it  was  decided  to  select  the  275-knot,  3- 

■ 

G  flight  position  as  the  maximum  loading  position  in  the 
normal  operating  envelope. 
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V.    FINITE  ELEMENT  ANALYSIS 

A.   THEORY 

The  method  of  finite  element  analysis  is  a  means  of 
simulating  the  structural  behavior  of  a  continuous  physical 
system  by  a  discretized  representation  of  that  system. 
Structures  are  represented  by  discrete  node  points  which  are 
connected  by  structural  elements.  The  nodes  form  a  grid  which 
details  the  general  shape  of  the  structure  while  the  elements, 
although  they  appear  as  only  lines,  are  mathematically  given 
the  physical  properties  of  the  portion  of  the  structure  which 
they  are  there  to  represent.  A  physical  structure  is  thus 
transformed  into  a  mathematical  representation  for  the  purpose 
of  analyzing  some  behavior  of  the  structure.  This  analysis 
may  be  in  the  area  of  dynamic  response,  heat  transfer,  or,  as 
is  the  case  here,  static  loading  response.  This  method  of 
analysis  is  widely  proven  to  be  highly  accurate  and  has  been 
used  in  many  engineering  fields,  including  aerospace, 
automotive,  civil  and  mechanical  applications.  With  the 
increased  capability,  speed  and  data  storage  capacity  of 
microcomputers,  this  analysis  technique  is  no  longer  limited 
to  mainframe  applications,  as  was  the  case  a  few  years  ago. 
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1 .   Accuracy 

The  accepted  rule  of  thumb  in  finite  element  modeling 
is  that  the  use  of  more  node  points  results  in  a  more  accurate 
solution.  Convergence  tables  have  been  developed  which  show 
that  the  use  of  fewer  nodes  increases  the  stiffness  of  the 
model  .  The  main  drawback  to  using  a  large  number  of  nodes  is 
that  it  greatly  increases  computation  time  and  requires  larger 
amounts  of  storage  space  than  a  model  of  the  same  structure 
using  fewer  nodes  [Ref.  14].  The  finite  element  analysis 
software  used  for  this  thesis  is  called  MSC/pal  2®  and  is  a 
product  of  the  MacNeal-Schwendler  Corporation,  the  company 
which  also  creates  the  highly  respected  NASTRAN®  finite 
element  application  for  VAX/VMS  work  stations  and  mainframe 
operations.  The  accuracy  of  MSC/pal  2®  has  been  tested  and 
documented  by  the  manufacturer  [Ref.  15].  In  addition,  the 
manufacturer  recommends  simple  hand  calculations  to  verify 
that  results  obtained  for  a  given  model  are  reasonable  (i.e. 
within  the  same  order  of  magnitude).  This  was  done  by 
calculating  simple  beam  bending  stress  with  constant  area 
cross  sections  approximating  the  rib  legs  of  the  model. 
Results  established  that  the  finite  element  model  produced  a 
solution  which  was  well  within  expected  norms. 
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2.   Equations 

The  number  of  equations  to  be  solved  in  the  finite 
element  analysis  is  equal  to  the  number  of  deqrees  of  freedom 
in  the  model.  Each  node  point  has  six  degrees  of  freedom: 
three  translational  (in  each  of  the  (x,y,z)  coordinate 
directions,  and  three  rotational  (about  each  of  the  coordinate 
axes).  Stiffness  equations  are  generated  for  the  stiffness  of 
each  connecting  element,  based  on  the  specified  material 
properties  (Young's  modulus,  shear  modulus,  mass  density, 
tensile  yield  stress  are  specified)  and  the  geometric 
configuration  of  the  element.  Elements  may  take  the  form  of 
beams,  triangles,  quadrilaterals  and  others.  These  nodal 
stiffnesses  are  combined  to  form  a  system  stiffness  matrix, 
[K],  of  size  (NxN)  where  N  is  the  number  of  equations. 
Degrees  of  freedom  may  be  eliminated  by  setting  them  to  zero 
in  the  model  definition  phase  (a  way  of  applying  boundary 
constraints)  or  fully  retained  as  was  done  in  this 
application.  (Here,  boundary  constraints  were  applied  in  the 
load  as  discussed  in  Chapter  IV.  This  application  allows 
variation  of  the  boundary  displacements  from  one  load  to  the 
next  and  allows  recovery  of  reaction  forces  at  the  constrained 
nodes. )  Once  the  stiffness  matrix  has  been  formed,  the  static 
analysis  may  be  performed  according  to  the  following  equation: 

[K]{U}  =  {F} 
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where , 


•  {F}  =  column  vector  of  applied  loads  (Nxl) 

•  {U}  =  resultant  column  vector  of  nodal  displacements 

(Nxl) 


Gaussian  elimination  is  employed  to  solve  the  matrix  equation. 
In  large  models,  as  is  the  case  here,  matrix  partitioning 
takes  place  prior  to  solution. 

Once  the  displacements  are  known,  stress-strain 
relationships  are  employed  to  compute  stress  values  throughout 
the  structure.  These  stresses  are  available  to  the  user  in 
the  form  of  major  and  minor  principal  axis  stresses,  Von  Mises 
stress  concentrations  and  maximum  shear  stresses.  In 
addition,  output  of  displacements  and  rotations  are 
accessible. 

B.   STRUCTURAL  REPRESENTATION 

1.   The  Leading  Edge  Structure 

The  leading  edge  segment  considered  for  this  analysis 
was  the  port  wing,  center  section,  located  between  the  number 
three  and  four  engine  nacelles.  Figure  17  shows  a  cutaway 
view  of  the  structure.  The  segment  is  composed  of  12  vertical 
ribs  supporting  a  double  (inner,  outer)  skin.  The  outer  skin 
is  .040  inches  thick  and  the  inner  is  a  stamped  corrugation  of 
.016  inches.  Assembly  of  the  outer  and  inner  skins  provides 
a  series  of  ducts  approximately  .25  inches  in  height  and  two 
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Figure  17.  P-3  leading  edge  structure  [Ref  16]. 

inches  wide  each.   These  ducts  run  longitudinally  along  the 
inner  surface  of  the  leading  edge  and  are  continuous  over  the 
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spanwise  length  of  the  leading  edge.  The  purpose  of  this 
series  of  ducts  is  to  provide  a  channel  through  which  bleed 
air  may  travel  to  heat  the  leading  edge.  The  tube  seen 
extending  through  the  structure  delivers  the  bleed  air  to  the 
ducts.  The  leading  edge  segment  is  secured  to  the  front  spar 
of  the  wing  by  a  full-length  piano  hinge  at  the  bottom  edge 
and  a  screwed-down  spar  cap  flange  at  the  top.  This 
arrangement  provides  access  to  the  area  for  maintenance 
functions. 

2.   The  Model 

The  basic  finite  element  model  used  in  this  analysis 
was  obtained  through  translation  of  a  NASTRAN®  finite  element 
model  using  a  function  in  the  MSC/pal  2®  application  called 
NASPAL®.  NASPAL®  reads  the  NASTRAN®  text  file  for  the  model 
and  rewrites  it  in  the  format  used  by  MSC/pal  2®.  The 
NASTRAN®  model  file  was  obtained  from  Aerostructures,  Inc. 
through  contact  with  NAVAIR's  AIR-530  office.  Because  the 
NASTRAN®  model  consisted  of  2749  nodes,  it  was  necessary  to 
reduce  the  model  size  to  meet  the  MSC/pal  2®  limitation  of 
2000  nodes.  This  reduction  was  done  by  entering  a  set  of 
geometric  coordinates  during  the  NASPAL®  translation  and 
instructing  the  translator  to  consider  only  the  outboard  64 
inches  of  the  model.  In  effect,  the  leading  edge  section  was 
severed  between  WS  247  and  WS  256,  or  between  the  third  and 
fourth  ribs  from  the  left  end  as  shown  in  Figure  17.   The 
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remaining  portion  of  the  model  was  translated  from  the  fully 

defined  model.  This  approach  was  considered  to  provide  a  more 

accurate  solution  than  increasing  the  spacing  between  nodes, 

keeping  the  highest  available 

level  of  detail  in  the  model 

(again,   more   nodes   mean 

better  accuracy  for  the  same 

structure).   In  doing  so,  it 

was  necessary  to  modify  the 

inboard   end   rib   of   the 

structure  (WS  256)  since  the 

end   ribs   were   constructed 

differently    than    the 

intermediate  ribs.   The  end 

ribs  are  closed  in  the  front 

and   have   single,   riveted 

flanges   on   their   interior 

(Figure  18,   bottom),  while 

the  intermediate  ribs  (Figure 

18,  top)  have  an  open  front 

to  allow  access  for  the  bleed 


air  "pump  cap"  assembly  which 
leads   to   the   double   skin 


described   earlier. 


The 


intermediate  ribs  also  have 


double   flanges   as   shown. 


Figure  18 


Intermediate  and 
end  rib  detail. 
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These  changes  to  the  model  were  made  in  the  rib  at  WS  25  6  (now 
the  inboard  end  rib  of  the  model)  by  adding  node  points  and 
quadrilaterals  to  close  the  front  and  by  rearranging  the  beam 
element  properties  within  the  MSC/pal®  2  model  definition  text 
file.  In  this  way,  the  model  was  altered  to  represent  a 
shortened  leading  edge  segment  with  properly  defined  ribs. 
The  model  was  constructed  as  an  all-aluminum  structure,  and  as 
stated  in  Chapter  IV,  the  loads  applied  were  generated  for 
each  specific  wing  station  and  scaled  according  to  the 
spanwise  distance  between  the  node  points.  In  this  way,  the 
model  detailed  here  received  a  scaled-down  load  for  its 
scaled-down  size.  The  upper  and  lower  surface  views  of  the 
model  are  presented  in  Figures  19  and  20.  The  end  views  of 
the  model  are  shown  in  Figure  21. 

The  lower  hinge  is  replicated  in  the  model  by  leaving 
the  Y-axis  rotation  unrestrained.  The  upper  flange  of  the 
model  is  secured  in  all  six  degrees  of  freedom.  Displacements 
for  front  spar  twist  are  incorporated  in  the  upper  and  lower 
constraints  as  detailed  in  the  previous  Chapter.  All  of  these 
boundary  conditions  are  input  through  the  "Displacements 
Applied"  command  section  of  the  load  file. 

Another  difference  between  this  model  and  the  original 
NASTRAN®  code  developed  for  NAVAIR  is  that  the  original  did 
not  incorporate  stiffness  generation  in  the  skin  of  the  model. 
The  skin  thickness  of  the  NASTRAN®  model  was  .056  inches, 
which  is  the  sum  of  the  outer  and  inner  skin  thicknesses  but 
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Figure  19.  Upper  surface  of  finite  element  model 


Figure  20.  Lower  surface  of  finite  element  model 
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without 


the 


stiffness 
generation  enabled 
in  the  model 

definition  file, 
the  skin  would  have 
no  stiffness.  In 
effect,  it  would 
act  as  a  non- 
loadbearing 
membrane.  In  the 
MSC/pal  2®  model, 
stiffness 
generation  was 

enabled,  but  when 
loaded,  this 

resulted  in 

excessive 
deformation   of    the 


skin. 


It       was 


decided   that   the 


.056   inch   thick 


Figure  21 


End  views  of  finite  element 
model  (outboard,  inboard). 


skin  did  not  accurately  model  the  combined  effect  of  the  inner 
and  outer  skin  combination  since  the  corrugated  inner  skin 
would  be  much  stiff er  than  its  mere  thickness  (.016  inches) 
would  represent.    Although  there  is  a  variable  stiffness 
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factor  in  the  MSC/pal  2  application,  available  references  on 
finite  element  analysis  did  not  detail  the  calculation  of  this 
factor.  Therefore,  it  was  decided  to  increase  the  skin 
thickness  to  an  equivalent  thickness  which  would  accurately 
replicate  the  behavior  of  the  corrugated  skin  combination. 
This  thickness  was  calculated  by  equating  the  moment  of 
inertia  of  a  box  beam  formed  by  a  single  corrugation  and  outer 
skin  combination  with  that  of  a  solid  cross  section  about  the 
same  reference  axis  (at  the  upper  surface).  After  allowing 
for  spaces  between  the  corrugations  where  the  two  skins  are 
riveted  together,  an  equivalent  skin  thickness  of  .1838  inches 
was  determined  and  employed.  This  action  reduced  the 
deformation  of  the  skin  under  load  to  reasonable  norms. 

Loads  were  applied  to  the  model  at  spanwise  rows  of 
selected  node  points  which  most  nearly  corresponded  to  the  mid 
points  of  the  panels  in  the  two-dimensional  panel  method. 
These  point  loads  were  applied  in  the  vertical  (Z)  and 
horizontal  (X)  directions  in  units  of  pounds  force. 

Material  properties  used  to  represent  the  2024-T6 
aluminum  structure  in  the  construction  of  the  model  were  as 
follows: 

•  Young's  modulus  (E)  =  1.06E+07  psi  [Ref.  17] 

•  Shear  modulus  (G)  =  4.0E+06  psi  [Ref.  17] 

•  Poisson's  ratio  (u)=  3.25E-01 

•  Tensile  yield  strength  (a  )  =  4.7E+04  psi  [Ref.  18] 
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where  Poisson's  ratio  (u)  was  calculated  according  to  the 
standard  formula: 

G  =  .  E 


2(l+v) 


C.   APPLICATION 

Many  load  conditions  were  examined  in  the  course  of  this 
thesis.  Presented  here  are  the  six  load  cases  which  give  the 
best  overall  illustration  of  the  observed  effects  of  static 
aerodynamic  loading,  both  with  and  without  static  aeroelastic 
effects  included.  In  this  way,  the  effect  of  wing  box  twist 
may  be  seen,  along  with  the  combined  effect  of  angle  of  attack 
and  dynamic  pressure.  In  Table  3,  the  features  of  the  six 
load  cases  are  given.  The  L/R  column  provides  a  distinction 
between  a  wings  level  (L)  pullup  or  a  rolling  (R)  pullup.  The 
effect  of  rolling  into  a  turn  during  the  application  of  G 
loading  (as  in  a  climbing  breakaway  maneuver)  is  not  the  same 
as  that  of  rolling  out  of  a  turn  during  G  application  (as  in 
rolling  to  wings  level  while  pulling  out  of  a  dive).  Since  it 
was  found  that  the  loading  effect  in  terms  of  both  twist  and 
air  loading  was  greater  during  the  former  (due  to  the  combined 
effect  of  roll  helix  angle,  aileron  deflection  and  air  load), 
the  former  was  chosen  for  presentation  here  in  the  rolling 
load  cases.  All  load  cases  were  generated  at  135,000  pounds 
gross  weight  except  for  number  4  which  was  done  at  110,000 
pounds.    In  load  case  number  1,  the  effect  of  twist  was 
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eliminated  from  the  load  solution  by  setting  the  MSC/pal  2 
displacements  to  zero  so  that  the  effect  of  wing  torsional 
twist  may  be  seen  by  comparison  with  load  case  2 . 

TABLE  3.   LOAD  CASES  EMPLOYED 

Load  Case  Airspeed  (kts)  Load  Factor  L  or  R      Comment 

1  275  3.0  Level     No  twist. 

2  275  3.0  Level 

3  275  2.4  Roll 

4  240  3.0  Level      110,000  # 

5  350  2.4  Roll 

6  325  4.5  Level 

1.   Input  Data 

Presented  in  Table  4  are  the  input  data  sets  for  each 
of  the  load  cases.  These  are  given  for  the  three  pertinent 
non-dimensional  spanwise  wing  stations  (2y/b)  as  generated  by 
the  span  load  program.  The  24  actual  wing  stations  (2y/b  = 
0.4322  through  0.5397)  employed  for  the  two-dimensional  panel 
program  were  solved  by  curve  fitting  the  lift  coefficients 
bounded  by  these  extremes,  as  described  in  Chapter  IV.  The 
92-inch  loads  are  the  approximate  longitudinal  (positive  aft) 
and  vertical  (positive  up)  total  loads  which  would  be  seen  by 
a  complete  center  section  leading  edge  segment  on  the 
aircraft.  The  total  load  applied  to  the  shortened  finite 
element  model  would  be  some  69.6  percent  of  this  stated  load. 
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TABLE  4.   FINITE  ELEMENT  ANALYSIS  INPUT  DATA 


Load  Case 

2y/b 

Cj0 

Twist  (deer) 

92"  Load  fibs) 

1 

.35 

1.2726 

0.0000 

FX 

= 

-5157 

.45 

1.2379 

0.0000 

FZ 

= 

13155 

.55 

1.1834 

0.0000 

2 

.35 

1.2726 

0.8736 

FX 

= 

-5157 

,45 

1.2379 

1.1221 

FZ 

= 

13155 

.55 

1.1834 

1.3764 

3 

.35 

1.0778 

0.7605 

FX 

= 

-3657 

.45 

1.0569 

0.9966 

FZ 

= 

11069 

.55 

1.0090 

1.2572 

4 

.35 

1.3465 

0.7216 

FX 

= 

-4454 

.45 

1.3118 

0.9272 

FZ 

= 

10666 

.55 

1.2575 

1.1375 

5 

.35 

0.7109 

0.5673 

FX 

= 

-2103 

.45 

0.6723 

0.7335 

FZ 

= 

10540 

.55 

0.6100 

0.9169 

6 

.35 

1.3789 

1.2997 

FX 

= 

-5982 

.45 

1.3247 

1.6613 

FZ 

= 

16742 

.55 

1.2453 

2.0279 

2.   Finite  Element  Analysis  Results 
a.      Load  Case   1 

In  looking  at  the  275  knot  load  without  static 
aeroelastic  twist,  it  was  noted  that  the  largest  stress 
concentrations  in  the  structure  were  located  in  the  rib  legs, 
with  the  lower  legs  experiencing  approximately  15  ksi  in 
tension  along  their  upper  flanges  and  the  upper  legs  seeing 
about  -14  ksi  (compression)  along  the  lower  flanges.  These 
stresses  were  evenly  distributed,  as  may  be  seen  by  the  values 
in  Figures  22,  23  and  24  where  the  major  principal  axis  stress 
(a-j-)  contours  are  shown.  The  minor  principal  (olz)  stresses, 
Von  Mises  and  shear  stresses  show  a  similarly  even 
distribution,  with  all  stress  levels  well  below  the  yield 
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stress  value  (o  s  =  47  ksi)  for  the  material.  It  should  be 
noted  that  the  apparent  deformations  in  the  plot  are 
exaggerated  and  not  to  scale.  This  scaling  provides  the 
viewer  with  a  better  perception  of  the  direction  of 
displacement  occurring,  although  in  reality,  the  displacements 
are  only  on  the  order  of  0.06  inches  for  this  load  case  as 
determined  by  MSC/pal  2®. 
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Figure  22.  Major  principal  axis  stress  contours  on 
inboard  (WS256)  end  rib  in  untwisted 
condition  (Case  1). 
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Major  principal  axis  stress  contours  on 
middle  (WS282)  rib  in  untwisted  condition 
(Case  1). 
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Figure   24.    Major  principal   axis   stress   contours   on 

outboard    (WS320)    rib   in   untwisted   condition 
(Case    1). 
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The  untwisted  load  case  (case  1)  does  not 
accurately  reflect  the  behavior  of  the  leading  edge  since  it 
does  not  include  the  effect  of  aeroelastic-induced  spanwise 
twisting  (torsion)  of  the  wing.  The  stress  contours  in 
Figures  22,  23  and  24  were  presented  for  comparison  to  load 
case  2. 

b.      Load  Case  2 

A  dramatic  difference  in  the  observed  stress 
contours  occurred  when  the  finite  element  model  was  subjected 
to  the  net  twist  occurring  in  the  wing  box  as  determined  by 
the  static  aeroelastic  span  load  analysis.  It  should  be  noted 
that  the  displacement  input  to  the  model  equated  to  only  about 
0.3  degrees  of  front  spar  rotation  from  WS  256  to  WS  320.  The 
same  set  of  stress  contour  plots  as  above  are  given  in  Figures 
25,  26  and  27.  A  striking  contrast  existed  between  the  first 
and  second  load  cases,  with  the  case  2  major  principal  axis 
stress  distribution  very  unevenly  spread  between  the  inboard 
and  outboard  ends.  As  may  be  seen  in  Figure  25,  the  inboard 
end  rib  (WS  256)  experienced  more  than  double  the  stress 
concentration  in  its  lower  leg  with  34.8  ksi  being  the  highest 
level  contour  shown.  Moving  outboard,  WS  282  (Figure  26) 
showed  a  maximum  of  about  20  ksi  in  its  lower  leg  while  the 
stress  in  the  same  leg  on  the  outboard  (WS  320)  rib  went  to 
zero.  (A  stress  level  of  2.9  ksi  in  the  saddle  and  upper  leg 
of  the  rib  may  still  be  seen. ) 
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Major  principal  axis  stress  contours  on 
inboard  (WS256)  end  rib  with  twist  applied 
(Case  2) . 
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Figure   26.    Major  principal *axis   stress   contours   on 
middle    (WS282)    rib  with   twist   applied 
(Case    2) . 
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Figure  27.  Major  principal  axis  stress  contours  on 
outboard  (WS320)  rib  with  twist  applied 
(Case  2) . 

Since  it  was  observed  that  this  same  pattern  of 
stress  distribution  occurred  for  the  minor  principal  axis,  Von 
Mises  and  maximum  shear  stresses,  only  the  inboard  end  (WS 
256)  contours  are  shown  in  Figures  28,  29  and  30.  This 
pattern  is  to  be  expected  since  they  are  all  geometrically 
related.  The  major  and  minor  principal  stresses  occur  on 
planes  on  which  there  is  no  shear  stress  and  are  oriented 
perpendicularly  to  each  other,  while  the  maximum  shear  stress 
occurs  on  planes  which  are  at  angles  of  45  degrees  to  the 
principal  planes.  Von  Mises  stresses  (av)  are  derived  from  a 
criterion  known  as  the  Maximum  Distortion  Energy  Criterion  and 
may  be  found  from  the  eguation, 

ov2  =  ox2   -  QjCn  +  aIZ2 
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which  is  based  on  the  determination  of  the  energy  associated 
with  changes  in  shape  of  a  given  material.  These  relation- 
ships are  valid  under  the  assumption  of  a  plane  stress 
condition  in  the  material.  This  means  that  the  metal  is  thin 
in  comparison  to  its  other  dimensions  so  that  the  stresses 
across  the  thickness  of  the  metal  may  be  considered  as 
negligible.  This  plane  stress  condition  is  the  case  for  most 
aircraft  structural  components  since  they  are  made  of  sheet 
material,  and  is  valid  here.  [Ref.  19,  20] 
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Figure  28.  Minor  principal  axis  stress  contours  on 

inboard  (WS256)  end  rib  with  twist  applied 
(Case  2) . 
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Figure   29.    Von  Mises   Criterion   stress   contours   on 

inboard    (WS256)    end   rib  with  twist   applied 
(Case   2) . 
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Figure   30.    Maximum  shear   stress   contours   on   inboard 

(WS256)    end   rib  with  twist   applied    (Case   2). 

In  Figures  31  through  36,    the  x  and  z  displacements  of 
the    ribs   are    shown    in   pairs    from   inboard   to   outboard.      Note 
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the  small  displacements  occurring  in  the  upper  flange  and 
lower  hinge  due  to  twist  of  the  front  spar. 
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Figure  31.  Displacement  (x)  at  inboard  rib  (WS256). 
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Figure    32.    Displacement    (z)    at   inboard   rib    (WS256). 


67 


ccc 


TRANS.   DEFL.   X 

R-3.6889E-8Z 
B-Z.7G67E-BZ 
C-1.8444E-BZ 
D-9.ZZZZE-83 
E  B.BB88E+88 
F  9.ZZZZE-B3 
G  1.8444E-8Z 
H  Z.YGG7E-BZ 
I  3.6889E-BZ 
J  4.6111E-BZ 


Figure   33.    Displacement    (x)    at  middle   rib    (WS282). 
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Figure   34.    Displacement    (z)    at  middle   rib    (WS282). 
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Figure  35.  Displacement  (x)  at  outboard  rib  (WS320). 
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Figure   36.    Displacement    (z)    at   outboard   rib    (WS320). 

Another  available  representation  of  the  stress 
concentrations  present  in  the  structure  is  an  X-Y  plot  of  all 
four  stresses  together  in  a  form  resembling  a  frequency 
scatter    on    an    oscilloscope.       This    plot    for    the    case    2    load 
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condition  is  presented  in  Figure  37.  While  this  is  a  rather 
cluttered  plot  it  does  serve  to  provide,  at  a  glance,  the 
maximum  and  minimum  stress  concentrations  in  the  structure. 
Since  the  stress  contours  of  individual  members  in  the 
structure  have  been  examined  and  it  has  been  found  that  the 
highest  stresses  are  in  the  rib  legs,  the  peak  stresses 
depicted  may  be  attributed  to  these  locations.  The  positive 
(upper)  portion  of  the  plot  is  a  combination  of  the  major 
principal,  Von  Mises  and  shear  stresses  while  the  lower  half 
depicts  the  minor  principal  axis  stresses.   For  the  sake  of 
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Figure  37.  X-Y  plot  of  combined  stresses  (Case  2). 

brevity,  and  since  the  contour  plots  of  the  individual  ribs 
retain  the  same  relative  form  as  those  presented  for  cases  1 
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and  2,  the  stress  levels  for  the  remaining  load  cases  will  be 
presented  in  this  form. 
c.  Load  Case  3 

For  the  275-knot,  2.4-G  rolling  pullup,  the 
maximum  stress  levels  observed  are  approximately  33  ksi  on  the 
descending  wing  as  shown  in  Figure  38.  The  loads  on  the 
leading  edge  of  the  ascending  wing  are  lower  than  those  for 
the  descending  wing.  This  speed  was  chosen  because  of  the 
relationship  which  was  shown  to  exist  between  leading  edge 
loading,  angle  of  attack  and  dynamic  pressure.  The  result  of 
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Figure  38.  X-Y  plot  of  combined  stresses  (Case  3). 
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higher  airspeed  at  the  same  aircraft  load  factor  may  be  see  in 
part  e  of  this  Chapter. 
d.      Load  Case   4 

The  110,000  pound  weight  condition  was  chosen  for 
this  load  case  in  order  to  examine  the  reduction  in  leading 
edge  stress  when  operating  within  the  normal  flight  envelope 
at  less  than  maximum  gross  weight.  The  240-knot  airspeed  is 
the  approximate  "corner"  speed  at  a  load  factor  of  3G  for  this 
weight.  The  maximum  stress  values  seen  under  this  condition 
are  28.8  ksi  in  tension  and  28.3  ksi  in  compression  as  seen  in 
Figure  39. 
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Figure  39.  X-Y  plot  of  combined  stresses  (Case  4). 
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e.  Load  Case   5 

In  load  case  5,  the  effect  of  a  2.4-G  rolling 
pullup  at  135,000  pounds  gross  weight  at  a  true  airspeed  of 
350  knots  was  examined.  In  Figure  40,  it  may  be  seen  that  the 
maximum  stress  values  are  approximately  26  ksi,  which  is 
approximately  7  ksi  less  than  the  same  maneuver  at  275  knots 
( load  case  3 ) . 
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Figure  40.  X-Y  plot  of  combined  stresses  (Case  5). 

f.      Load  Case   6 

In  order  to  examine  the  design  ultimate  flight 
load  condition,  load  case  6  was  run  at  325  knots  and  4.5G. 
The  airspeed  corresponds  to  the  approximate  stall  speed 
(corner  speed)  for  the  135,000  pound  airplane  gross  weight. 
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Stress  values  of  nearly  52  ksi  may  be  seen  in  Figure  41.  Note 
that  this  result  exceeds  the  yield  stress  for  the  material  (47 
ksi)  . 
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Figure  41.  X-Y  plot  of  combined  stresses  (Case  6). 

3.   Discussion 

The  primary  factor  of  concern  in  the  observed  loads 
presented  above  is  the  apparently  high  level  of  stress  in  the 
rib  legs  as  determined  by  the  methods  stated  in  the  course  of 
this  writing.  According  to  the  Engineering  Investigation  (EI) 
report  for  the  USN  P-3  mishap  [Ref.  21]  and  a  preliminary 
report  issued  by  Aircraft  Research  Laboratories  of  Melbourne, 
Australia  [Ref.  22],  the  initial  structural  failure  in  both 
mishaps  was  believed  to  have  occurred  in  the  lower  rib  legs, 
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with  initial  failure  occurring  in  the  outboard  leg  and 
proceeding  inboard  in  a  series  of  sequential  failures.  It 
would  appear  from  the  results  of  this  study  that  this  failure 
sequence  may  have  been  reversed,  since  the  highest  observed 
stresses  were  in  the  inboard  rib.  Additionally,  the  fact  that 
there  is  an  intermediate  rib  at  WS  317,  only  three  inches  from 
the  outboard  end  rib  (rather  than  the  eight-  to  nine-inch 
spacing  for  all  other  ribs  in  the  leading  edge)  supports  the 
idea  that  the  failure  would  likely  have  initiated  at  some 
other  location.  This  close  proximity  of  two  ribs  means  that 
the  load  in  that  region  would  be  shared,  thereby  reducing  the 
stress  in  the  outboard  rib.  In  any  case,  the  primary  question 
remains  as  to  whether  there  were  sufficiently  high  stress 
levels  in  the  rib  legs  to  cause  structural  failure. 

While  the  maximum  observed  stress  levels  (found  here 
for  operation  within  the  flight  envelope,  occurring  under  load 
case  2  in  a  wings  level,  3-G  pullup  at  135,000  pounds  gross 
weight)  were  some  12  ksi  below  the  yield  stress  for  the 
material,  two  areas  of  concern  must  be  addressed:  1)  the 
additional  effect  of  extending  this  data  from  the  assumed  64- 
inch  model  size  to  the  true  92-inch  leading  edge  segment  which 
is  installed  on  the  aircraft;  and  2)  the  effect  of  stress 
concentrations  around  holes  in  the  rib  legs. 
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a.  Extending  the  Structure 

Proceeding  under  the  assumption  that  the  stress 
concentrations  vary  linearly  with  the  length  of  the  structure, 
extending  the  64-inch  model  to  the  full  (92-inch)  length  of 
the  true  structure  would  result  in  an  increase  in  the  observed 
stress  in  the  lower  rib  leg  to  51.9  ksi,  which  is  in  excess  of 
the  47  ksi  yield  stress.  Examining  the  major  principal  axis 
stresses  in  the  lower  legs  of  WS  256  and  WS  282,  this 
linearity  would  appear  to  hold  true.  Since  aluminum  is  a 
ductile  material,  the  stresses  may  subside  as  the  material 
begins  to  yield,  therefore,  the  structure  may  still  carry  the 
load  without  fracture.  However,  changes  in  the  pressure 
distribution  around  the  leading  edge  must  also  be  considered 
as  the  structure  begins  to  deform.  Given  the  observed 
directions  of  the  deformations  as  seen  in  the  contour  plots, 
it  seems  likely  that  the  aerodynamic  load  would  continue  to 
increase  due  to  the  increased  camber  in  the  leading  edge  as  it 
lifts  perpendicular  to  the  chord  and  develops  a  more  circular 
form.  This  deformation,  although  initially  small  could  add  to 
the  load,  initiating  a  divergent  stress  scenario  leading  to 
failure. 

b.  Stress  Concentrations 

It  is  an  established  fact  that  stress 
concentrations  around  circular  holes  can  multiply  the 
localized  stress  by  factors  of  between  two  and  four,  depending 
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on  the  directional  orientation  of  the  loading  seen  by  the 
specimen.  Along  the  inner  surface  of  the  end  ribs  is  a 
continuous  flange  which  is  attached  by  a  rivet  spacing  of 
approximately  one-half  inch.  Additionally,  there  are  four 
tapered  angle  stiffeners  riveted  to  the  web  section  of  the 
rib.  Two  of  these  are  symmetrically  located  on  the  upper  and 
lower  legs  at  a  position  approximately  52  percent  of  the 
leading  edge  chordwise  dimension  forward  of  the  hinge.  (On 
the  outboard  end  rib  (WS  320)  this  is  approximately  12  inches 
forward  of  the  hinge. )  Around  these  rivet  holes,  stress 
concentrations  may  be  assumed  to  occur.  The  occurrence  of 
these  concentrations  at  the  onset  of  stress  loading  depends  on 
the  rivet  installation  procedure;  i.e.,  the  stress 
concentration  effect  is  delayed  if  the  rivet  hole  is  placed  in 
compression  (as  in  a  dimpled  or  double  dimpled  installation 
[Ref.  18]).  If  the  hole  were  initially  in  compression,  a 
margin  or  stress  buffer  would  be  provided,  as  the  material  is 
subjected  to  tension  stress;  i.e.,  the  tensile  stress  applied 
must  first  exceed  the  compressive  pre-stress  of  the  rivet 
installation  before  the  hole  will  begin  to  experience  a  stress 
concentration.  The  type  of  installation  of  the  rivets  under 
consideration  is  unknown  by  this  author.  Since  the  inside 
flange  of  the  upper  legs  is  placed  under  compressive  loading 
nearly  equivalent  to  that  experienced  by  the  lower  leg  in 
tension,  the  effect  of  compressive  concentrations  must  also  be 
considered.     If   the   rivet   holes  were   pre-stressed   in 
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compression,  this  would  add  to  the  compressive  stress  level  at 
this  location,  potentially  causing  the  upper  leg  to  fail  first 
in  compression.  As  before,  this  all  depends  on  the  stress 
reduction  due  to  yielding  which  takes  place  in  the  ductile 
material. 

In  the  Australian  incident  report,  the  failure 
sequence  was  described  as  having  initiated  along  the  line  of 
the  riveted  stiffener  on  the  lower  leg  of  the  outboard  end  rib 
as  mentioned  above.  In  addition,  it  goes  on  to  say  that  the 
failure  then  proceeded  along  a  line  of  dimples  in  the  lower 
legs  of  the  intermediate  ribs,  from  outboard  to  inboard. 
These  dimples  act  as  stiffeners  for  the  intermediate  ribs  and 
would  also  act  as  stress  concentrations  under  load.  None  of 
these  stiffeners,  on  either  outboard  or  intermediate  ribs,  nor 
their  rivet  holes,  are  included  in  the  finite  element  model 
employed  in  this  study.  Again,  the  failure  sequence  described 
is  opposite  of  that  expected  from  the  results  of  this  study. 
The  RAAF  determination  of  the  failure  sequence  was  based  on 
the  observed  outboard-to-inboard  bending  of  the  rib  fragments 
which  remained  attached  to  the  structure.  [Ref.  22] 
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VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.   CONCLUSIONS 

The  results  of  this  study  indicate  that  failure  of  a  P-3 
wing  leading  edge  segment  could  be  predicted  to  occur  within 
the  normal  aircraft  operating  envelope  providing  proper 
account  is  given  to  the  influences  of  wing  static  aeroelastic 
effects  both  upon  wing  span  loads  and  torsional  twists  induced 
in  the  wing  spar  box.  The  location  of  the  highest  stresses 
are  in  the  area  of  the  observed  failures  as  reported  in  the  EI 
and  Preliminary  reports  [Ref.  21  and  22].  It  is  reasonable  to 
assume  that  stress  concentrations  around  the  rivet  holes  could 
be  on  the  order  of  1.5  or  2.0  times  the  observed  stresses. 
This,  combined  with  the  effect  of  extending  the  model  to  the 
full  92-inch  length,  could  result  in  stress  levels  in  excess 
of  the  ultimate  strength  of  the  material,  60  ksi  [Ref.  18]. 

The  primary  evidentiary  conflict  with  this  study  is  the 
reported  location  of  apparent  failure  initiation  at  the 
outboard  vice  inboard  end.  The  observed  inboard  bending  of 
the  rib  fragments  may  have  been  induced  by  the  upper  portion 
departing  the  wing  in  some  fashion  other  than  that  assumed. 
It  is  difficult,  at  best,  to  assess  the  direction  in  which  the 
ribs  would  bend  as  the  failure  progressed. 
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While  fatigue  is  potentially  a  contributor  to  this  mode  of 
failure,  given  the  possibility  of  plastic  strain  having 
occurred  earlier  under  a  reduced  load  condition,  the  EI  [Ref. 
21]  stated  that  no  undue  hardness  was  detected  in  the  failed 
USN  structure.  This  would  indicate  that  no  material  strain 
hardening  had  taken  place  up  to  the  time  of  failure. 

Corrosion  is  another  potential  problem,  particularly  when 
dealing  with  aircraft  operating  routinely  in  low  level 
overwater  environments;  however,  no  levels  of  corrosion  which 
would  effect  the  structural  integrity  of  the  leading  edge  were 
cited  in  the  EI. 

Whether  or  not  the  aircrews  overstressed  the  aircraft  by 
exceeding  the  limitations  of  the  flight  envelope  cannot  be 
concluded  here.  Even  they  may  not  know  for  sure  whether  this 
was  the  case.  Flight  experiments  conducted  in  the  2F87F  P-3 
simulator  at  NAS  Moffett  Field,  CA,  during  the  course  of  this 
analysis  indicate  that  at  high  gross  weights  and  aft  center  of 
gravity  conditions,  it  is  easy  to  exceed  the  3G  limit  with 
just  a  firm  pull  on  the  control  yoke.  The  rapid  onset  of  G 
overload  may  be  imperceptible  due  to  the  rate  at  which  it  was 
observed  to  occur.  The  cockpit  G  meter  is  not  an  instrument 
which  is  normally  kept  in  the  pilot's  scan. 
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B.   RECOMMENDATIONS 

This  analysis  should  be  confirmed  by  independent 
validation.  If  these  results  are  verified,  re-enforcement  of 
the  end  rib  may  be  feasible;  alternately,  some  limitations 
should  be  placed  on  the  operational  flight  envelope  of  the 
P-3 .  Covered  here  are  some  recommended  procedures  for 
validation  along  with  approximations  of  interim  flight 
envelope  limits,  should  validation  and/or  repair  not  be 
possible  in  the  near-term. 

1.   Validation 

Since  the  Naval  Air  Systems  Command,  through 
Aerostructures,  Inc.,  has  the  base  model  from  which  this 
finite  element  model  was  derived,  it  is  recommended  that 
validation  be  conducted  with  the  following  modifications  on 
the  existing  model: 


•  Skin  stiffness  should  be  included  in  the  model.  The 
contribution  of  the  skin  to  the  load  bearing  properties  of 
the  structure  is  essential.  Skin  stiffness  could  be 
simulated  by  artificially  increasing  the  thickness  to 
simulate  the  double  skin  as  was  done  here,  or  preferably, 
through  addition  of  a  w  axis  (quadrilateral  element  local 
coordinate  system)  stiffness  factor  which  would  emulate 
the  combined  stiffness  of  the  inner  and  outer  skins. 

•  Static  aeroelastic  twist  must  be  included  in  the  analysis 
since  it  was  the  factor  contributing  to  the  dramatic 
increase  in  inboard  end  rib  stress.  Since  the  NASTRAN® 
model  is  already  of  the  appropriate  length,  this  model 
would  serve  to  verify  or  discount  the  theory  that  the 
stress  concentrations  would  increase  linearly  with 
extended  length. 

•  A  small  finite  element  model  of  the  area  of  immediate 
concern   (the   lower   leg  of   an   end   rib)   should   be 
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constructed  and  analyzed.  This  model  should  be  of 
sufficient  detail  to  include  the  rivet  holes  and 
stiffeners  in  order  to  examine  the  amount  of  stress 
concentration  occurring.  Any  pre-stress  of  the  rivet 
holes  should  be  taken  into  account. 


Laboratory  tests  of  actual  end  ribs  would  prove 
helpful  in  determining  the  effect  of  the  rivet  holes  on  stress 
concentrations  and  provide  data  on  the  load  levels  required  to 
cause  fracture. 

2.   Flight  Envelope  Modifications 

If  these  results  are  verified,  the  P-3  flight  envelope 
should  be  restricted  to  prevent  another  mishap.  If,  in  the 
estimation  of  NAVAIR,  verification  cannot  be  accomplished 
within  a  reasonable  amount  of  time  (six  months),  then  it  is 
recommended  that  an  interim  measure  be  taken  to  restrict  the 
flight  envelope  until  such  verification  or  disproval  can  be 
accomplished.  While  the  precise  envelope  restrictions  would 
be  developed  by  NAVAIR,  it  is  recommended  on  the  basis  of  the 
results  observed  here  that  the  following  limitations  be 
applied  if  interim  limitations  are  deemed  appropriate: 


Reduce  the  maximum  sustained  load  factor  for  the  aircraft 
from  3G  to  2G  at  all  operating  weights  above  100,000 
pounds . 

Reduce  the  maximum  sustained  load  factor  for  the  aircraft 
from  3G  to  2.5G  at  all  operating  weights  up  to  and 
including  100,000  pounds. 
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While  these  are  rough  estimates,  they  are  believed  to  provide 
a  sufficient  margin  of  safety  to  allow  safe  flight  without 
severely  impacting  daily  P-3  operations. 
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APPENDIX  A 


10   '  Program  "P3SPNLD3 .BAS" 

Solve  Static  Aeroelastic  Spanload  Problem 

for  a  wing  with  straight  taper 
Vortex  Lattice  capability  uses  swept  bound  vortex 

elements 
Inputs  consist  of: 

AR  =  Wing  Aspect  Ratio  =  B"2/S 

TR  =  Wing  Taper  Ratio  =  Ct/Cr 

SWP25  =  0.25  ch'd.  sweep  angle,  +'ve  is  sweepback 

CEA  =  Elastic  axis  location  on  const,  fraction 

chord  line 
MACH  =  Subsonic  Mach  no.  for  aerodynamic  compress. 

correctn. 
M   =  No.  of  equal  length  spanwise  stations,  RH  wing 
N   =  No.  of  equal  chordwise  stations 
**  Comments  ** 
Allows  Vortex  Lattice  solution  with  MxN  boxes  on 

RH  wing 
Max.  of  MxN  =50:'  Consistent  with  dimension 

statements 
N  =  1  for  elementary  (Modified  Weissinger)  lifting 

line  theory 
Wingspan,  B,  set  to  1188  inches  (for  P3  applicatn. ) 
100  '  $DYNAMIC 

DIMX1(50),  Yl(50),  X2(50),  Y2(50),  X3(50),  Y3(50), 
SWP(50) ,DIM  Al(50,  50),  A2(50,  50),  ASYM(50,  50), 
AANT(50,  50), DIM  S(50,  50),  F(20,  50),  G(20,  50), 
XEA(20),  YEA(20),  EI(20),  GJ(20),DIM  SWPM(50), 
DELA(50,  20),  MX(20,  50),  MY(20,  50 ) ,DIM  ALPHA( 50 ) , 
U(20,  20),  UGJ(20,  20),  UEI(20,  20),DIMA(50,  50), 
XIN(50):  '  for  Spanload  solutns.  using  ELU/SLVB  S/R's  . 
DIM  ALPHA1(50) ,  ALPHA2(50),  ALPHA3(50),  ALPHA4(50), 
Q(50),  DELCAM(5) 

OPEN  "C: QBFILES\EXTRA.DAT"  FOR  OUTPUT  AS  #1 
150  '  Input  wing  geometric  information  . .  P3  usage 
AR  =  7.539:  TR  =  .40088:  M  =  10 
'PRINT  :  PRINT  "Aspect  Ratio,  AR  =" ;  AR 
'PRINT  "Taper  Ratio,  TR  =";  TR 
'PRINT  "No.  R.H.  Wing  Spanwise  Stas.  =";  M 
'MACH  =  .6:  PRINT  TAB (5);  "Mach  No.  =";  MACH 
INPUT  "Airspeed  (knots)  =";  VEL 
MACH  =  (VEL  *  1.68781)  /  1116.3 
PRINT  "Mach  No.  =",  MACH 
IF  MACH  >  .95  THEN  GOTO  151 
FC  =  1!  /  SQR(1!  -  MACH  "2):  GOTO  152 
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151 

FC  =  1! 

152 

'  Continue 

N  =  5 

'PRINT  "No 

'INPUT  "No 

SWP25  =  -1 

Dummy  statement 

Chordwise  Stas.  =";  N 

Chordwise  Stas.  =";  N 
312:  'PRINT  "C/4  Sweep  (deg. )  =";  SWP25 
'INPUT  "0.25  chd.  sweep  (deg.)  =";  SWP25 
CEA  =  .4:  'PRINT  "Elast.  Axis  Locatn. ,  X/C  =";  CEA 
'INPUT  "Elastic  Axis  Locatn.,  X/C  =";  CEA 
'INPUT  "Dynamic  Pressure,  Q  (psi)  =";  Q 
INPUT  "Is  this  a  RIGID  run?  l=Yes,  2=No,  ANS=",  ANS 
IF  (ANS  =  1)  THEN 
Q  =  0 
ELSE 

Q  =  (.7  *  2116.2  *  MACH  A  2)  /  144 
END  IF 

PRINT  "Q  =",  Q 
KMAX  =  M  *  N 

'  Calculate  propeller  thrust  per  engine  (using  curve  fit 

'  from  2F87F  P-3  simulator  data) 

THRUST  =  (5849  -  5.069  *  VEL) 

'  Calculate  velocity  boost  due  to  prop  slipstream 

v  =  VEL  *  1.68781 

VI  =  (-V  +  SQR((v  A  2  +  (4  *  .05  *  THRUST  /  (143.14  * 

.0023769)))))  *  10 
'  Use  velocity  boost  to  find  boosted  Q  (Ql) 
Ql  =  (.5  *  .0023769  *  (v  +  VI)  ~  2)  /  144 
IF  (ANS  =  1)  THEN 

Ql  =  0 
END  IF 

'  Set  up  Q  as  a  vector  with  boosted  pressure  (Ql)  in 

'  region  behind  propellers 

FOR  I  =  1  TO  KMAX 

Q(I)  =  Q:  NEXT  I 

FOR  I  =  11  TO  35 

Q(I)  =  Ql:  NEXT  I 

PI  =  4!  *  ATN(1!):  '  Establish  constant 

SWP25  =  SWP25  *  PI  /  180!:  '  Convert  sweep  to  radians 

180  '  Print  header 
'GOSUB  1100 
'FOR  MM  =  1  TO  4 


200  '  Determine  Initial  Wing  Geometry 

B  =  1188!:    '  DEFAULT  Value  of  P3  Wing  Span,  inches 
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DELB  =  .5  *  B  /  M:  'Spanwise  spacing  increments  of 

'  vortex  elements 

CR  =  2!  *  B  /  (AR  *  (1!  +  TR) ) :  '  Root  Chord,  inches 

S  =  B  A  2  /  AR:  '  Wing  Area,  sg.  in. 

TANLE  =  TAN(SWP25)  +  ( . 5  *  CR  *  ( 1 !  -  TR) )  /  B: 

'  Tangent  L.E.  sweep 

'  Develop  Mean  Aero.  Chord  information 

MAC  =  2!  *  CR  *  ((1!  +  TR)  -  (TR  /  (1!  +  TR) ) )  /    31: 

'    Cmac ,  inches 

IF  TR  =  11  THEN  GOTO  210 

YMAC  =  .5  *  B  *  (1!  -  (MAC  /  CR) )  /  (1!  -  TR) :  GOTO  211 

210  YMAC  =  .5  *  B 

211  XMAC25  =  (YMAC  *  TANLE)  +  .25  *  MAC 


250  '  Determine  coords,  for 
'  control  pts. 
CONST1  =1!  -  TR 
FOR  I  =  1  TO  M 
CI  =  CR  *  (1!  -  (CONST1 
'  chord,  inbd 
C2  =  CR  *  (1! 
'  vortex 
C3  =  CR  *  (1! 
'  control  point 
DELC1  =  CI  /  N: 


vortex 

-  (CONST1 

-  (CONST1 
sta. 
DELC2  = 


wing  vortex  lattice  corners  & 

*  (I  -  1!)  /  M)):  'Total  wing 

*  I  /  M)):   '  Wing  chord,  outbd 
*(I-.5)/M)):  '  Chord  at 


C2  / 


YEA(I)  -  (I  - 
'  coordinates 
XEA(I)  =  YEA(I) 
FOR  J  =  1  TO  N 
K  =  (I  -  1)  *  N 
'  scheme 


5)  *  DELB: 


N:  DELC3  = 
Create  "M" 


C3  /  N 
Elastic 


axis 


*  TANLE  +  CEA  *  C3 


+  J:  'Create  Vortex 


( 


DELB:  Y2(K)  =  Y 

5  *  DELB) 
(J  -  .75)  +  Y1(K) 
(J  -  .75)  +  Y2(K) 
(J  -  .25)  +  Y3(K) 
-  X1(K))  /  (Y2(K) 


Y1(K)  =  (I  -  1) 

Y3(K)  =  Y1(K)  + 

X1(K)  -  DELC1  * 

X2(K)  =  DELC2  * 

X3(K)  =  DELC3  * 

TANSWP  =  (X2(K) 

SWP(K)  =  ATN( TANSWP) 

TANSWPM  =  FC  *  TANSWP 

SWPM(K)  =  ATN( TANSWPM) 

NEXT  J 

NEXT  I 

'  Tangent  of  Elastic  Axis 

DELEA  =  (XEA(2)  -  XEA(l)) 


X 

Lattice 

numbering 

(K)  +  DELE 

t 

* 

TANLE 

* 

TANLE 

* 

TANLE 
Y1(K)) 

Sweep,  etc. 
REA  =  SQR( DELEA 


2  +  DELB~2) 


TANEA  =  DELEA  /  DELB:  SINEA  =  DELEA  /  REA 

COSEA=DELB/REA 

WRITE  #1,  Q 

INPUT  "ENTER  G  LOAD,  n  =" ,  NZ 

'Calculate  required  CL 

CLREQD  =  (NZ  *  135000)  /  (.7  *  2116.2  *  MACH 

PRINT  "CLREQD  =",  CLREQD 


1300) 
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WRITE  #1,  NZ 
270  '  Bring  in  EI(I)  and  GJ(I)  values  for  Structural  model 
'1=1  (Root  Sta.)  &  I  =  M  (Tip  Sta.),  M  =  10  by 
'  default 
'  EI  Data,  Est.  for  P3  wing 


300 


DATA  8.30E+10,  6.80E+10,  5.27E+10,  4.10E+10,  2. 

,80E+10 

DATA  2.00E+10,  1.33E+10,  0.85E+10,  0.55E+10,  0. 

.35E+10 

FOR  I  =  1  TO  M:  READ  EI(I):  NEXT  I 

'  GJ  Data,  Est.  for  P3  wing 

DATA  7.50E+10,  5.25E+10,  3.90E+10,  2.90E+10,  2. 

,  10E+10 

DATA  1.30E+10,  0.85E+10,  0.65E+10,  0.35E+10,  0. 

•25E+10 

FOR  I  =  1  TO  M :  READ  G J ( I ) :  NEXT  I 

'Develop  Aerodynamic  Influence  Coefficients 

'  A1(K1,K2)  =  RH  Wing 

'  A2(K1,K2)  =  LH  Wing 

FOR  Kl  =  1  TO  KMAX:  '  Kl  Sta.  is  at  the  Control  Point 

FOR  K2  =  1  TO  KMAX:  '  K2  Sta.  is  at  the  Vortex 

Station 

NUM1  =  X3(K1)  -  XI (K2):  NUM2  =  X3(K1)  -  X2(K2) 

NUM1  =  FC  *  NUM1:  NUM2  =  FC  *  NUM2 

'  Bring  in  Prandtl-Glauert  factor 

DEN1  =  Y3(K1)  -  Y1(K2):  DEN2  =  Y3(K1)  -  Y2(K2) 

Rl  =  SQR(NUM1  A  2  +  DEN1  A  2 )  :  R2  =  SQR(NUM2"2+DEN2/S2 ) 

'  Find  trig,  functns  for  orthogonal  transformation  on 

'  swept  bound  vortex 

SINSWP  =  SIN(SWPM(K2) ) :  COSSWP  =  COS ( SWPM(K2 ) ) 

H  =  NUM1  *  COSSWP  -  DEN1  *  SINSWP 

Y1ROT  =  NUM1  *  SINSWP  +  DEN1  *  COSSWP 

Y2ROT  =  NUM2  *  SINSWP  +  DEN2  *  COSSWP 

COSTHET1  =  Y1ROT  /  Rl:  COSTHET2  =  -Y2ROT  /  R2 

'  Logic  Check  to  avoid  division  by  zero 

IF  (ABS(H))  <=  .001  THEN  GOTO  310 

DELWBD  =  (COSTHET1  +  COSTHET2 )  /  H:  GOTO  311 

310  DELWBD  =0!:    'No  downwash  contributn.  from  bound 
'  vortex 

311  '  Dummy  statement  space 

COS1  =  NUM1  /  Rl:  COS2  =  NUM2  /  R2 

DELWLH  =  (1!  +  COS1)  /  DEN1 

DELWRH  =  (1!  +  COS2)  /  DEN2 

A1(K1,  K2)  =  (DELWLH  +  DELWBD  -  DELWRH)  /  (8!  *  PI) 

NEXT  K2 

NEXT  Kl 

'  **  Similar  logic  for  LH  Wing  Panel  Aero.  Influence 
'     coeffs. 

315  FOR  Kl  =  1  TO  KMAX:  '  Kl  Sta.  is  at  the  Control  Point 
FOR  K2  =  1  TO  KMAX:  '  K2  Sta.  is  at  the  Vortex  Station 
NUM2  =  X3(K1)  -  X1(K2):  NUM1  =  X3(K1)  -  X2(K2) 
NUM2  =  FC  *  NUM2:  NUM1  =  FC  *  NUM1 
'  Bring  in  Prandtl-Glauert  factor 
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DEN2  =  Y3(K1)  +  Y1(K2):  DEN1  =  Y3(K1)  +  Y2(K2) 

Rl  =  SQR(NUM1  A  2  +  DEN1  ~  2):  R2  =  SQR(NUM2~2+DEN2~2 ) 

'  Find  trig,  functns  for  orthogonal  transformation  on 

'  swept  bound  vortex 

SINSWP  =  -SIN(SWPM(K2) ) :  COSSWP  =  COS ( SWPM(K2 ) ) 

H  =  NUM1  *  COSSWP  -  DEN1  *  SINSWP 

Y1ROT  =  NUM1  *  SINSWP  +  DEN1  *  COSSWP 

Y2ROT  =  NUM2  *  SINSWP  +  DEN2  *  COSSWP 

COSTHET1  =  Y1ROT  /  Rl:  C0STHET2  =  -Y2R0T  /  R2 

'  Logic  Check  to  avoid  division  by  zero 

IF  (ABS(H))  <=  .001  THEN  GOTO  320 

DELWBD  =  (COSTHET1  +  COSTHET2 )  /  H:  GOTO  321 

320  DELWBD  =0!:    'No  downwash  contributn.  from  bound 
'  vortex 

321  '  Dummy  statement  space 

COS1  =  NUM1  /  Rl:  COS2  =  NUM2  /  R2 

DELWLH  =  (1!  +  COS1)  /  DEN1 

DELWRH  =  (1!  +  COS2)  /  DEN2 

A2(K1,  K2)  =  (DELWLH  +  DELWBD  -  DELWRH)  /  (8!  *  PI) 

ASYM(K1/  K2)  =  A1(K1,  K2 )  +  A2(K1,  K2 ) 

AANT(K1,  K2)  =  A1(K1,  K2 )  -  A2(K1,  K2 ) 

NEXT  K2 

NEXT  Kl 

400  '  Determine  Struct.  Infl.  Matrix,  based  on  test  of  Q  >  0 
'  Note:  Q  =  0  psi  case  is  rigid  wing  situation 
IF  Q  >  0!  THEN  GOSUB  3000 

'  Skip  the  option  for  Sym.  or  Anti-Symm.  Soln. 
'GOTO  800:  'Symm.  solution  branch 

700  '  Select  Symmetric  or  Antisymmetric  Solution 

INPUT  "Symm.  or  AntiSymm.  Prob.,  S/A";  P$ 
710  IF  P$  =  "S"  THEN  GOTO  800 

IF  P$  =  "A"  THEN  GOTO  900 

GOTO  1000 
800  '  Find  Additional  Loading 

'  GOTO  805:  '  Branch  to  find  effect  of  wing  washout 

'  Introduce  Additional  type  of  alpha 

INPUT  "ENTER  AOA  IN  RADIANS,  AOA  =  ",  AOA 

FOR  K  =  1  TO  KMAX:  ALPHA(K)  =  AOA:  NEXT  K 

'GOTO  805 

'PRINT  "Elastic  wing,  Alpha  =  ",  ALPHA(l) 

WRITE  #1,  ALPHA(l) 

'  Introduce  alpha  due  to  built-in  geometric  twist 
'  0.0  deg.  at  root,  linearly  to  -2.5  deg.  at  tip 
WASH  =  -2.5  *  PI  /  180! 
FOR  I  =  1  TO  M 

ETA  =  (I  -  .5)  /  M:  DELWASH  =  ETA  *  WASH 
FOR  J  =  1  TO  N 
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K=(I-1)*N+J:  '  Create  panel  numbering 
ALPHAl(K)  =  DELWASH:  NEXT  J 
NEXT  I 
'GOTO  805 

'  Introduce  alpha  (rad.)  due  to  eng./ProP  assbly. 

'  dead-wgt . 

DATA  -0.0012,  -0.0042,  -0.0084,  -0.0138,  -0.0189 

DATA  -0.0240,  -0.0318,  -0.0366,  -0.0366,  -0.0366 

'  Read  dead-wgt .  alphas  one  at  a  time . . 

FOR  I  =  1  TO  M 

READ  DELALPH 

FOR  J  =  1  TO  N 

K=(I-1)*N+J:  '  Create  panel  numbering 

ALPHA2(K)  =  NZ  *  DELALPH:  NEXT  J 

NEXT  I 
t 

'    Introduce  alpha  (rad.)  due  to  camber 

DATA  -.0525,  -.0240,  0.0,  .0420,  .0730 

'  Read  in  alphas  due  to  camber  at  the  root  (WS  =  0) 

FOR  I  =  1  TO  N 

READ  DELCAM(I) 

NEXT  I 

'  Calculate  the  fifty  slopes  (one  per  control  point)  as 

'  airfoil  varies  from  14%  to  12%  thickness  from  root  to 

'  tip 

FOR  I  =  1  TO  M 

FOR  J  =  1  TO  N 

K=(I-1)*N+J 

ALPHA3(K)  =  DELCAM(J)  *  ( 1!  -  (I  -  .5)  /  (7  *  M) ) 

NEXT  J 

NEXT  I 

'  Introduce  alpha  (rad.)  due  to  aileron  float  angle 

'  Calculate  aileron  float  angle  (rad)  based  on  velocity 

'  (lets) 

FLAIL  =  -(1.4548  -  .0090025  *  VEL  +  .0000444  *  VEL  A  2) 

*  .01745 
FOR  I  =  1  TO  KMAX 
ALPHA4 ( K )  =  0 ! :  NEXT  I 
ALPHA4(39)  =  .2  *  FLAIL 

ALPHA4(44)  =  .2  *  FLAIL:  ALPHA4(49)  =  .2  *  FLAIL 
ALPHA4(35)  =  .35  *  FLAIL 
ALPHA4(40)  =  FLAIL:  ALPHA4(45)  =  FLAIL:  ALPHA4 ( 50 )=FLAIL 

805  '  Combine  angles  by  choice  for  net  alpha  distrib. 

INPUT  "Enter  type  of  run  desired:  0  =  Additional  only, 
1=  Washout  only,  2  =  Dead  wt  only,  3  =  Camber  only, 
4=  Aileron  float  only,  5  =  Complete  solution,  TYPE  =",  T 
WRITE  #1,  T 
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FOR  I  =  1  TO  KMAX 

IF  (T  =  0)  THEN 

ALPHA(I)  =  ALPHA(I):  'Additional  loading  selection 

ELSEIF  (T  =  1)  THEN 

ALPHA(I)  =  ALPHAl(I):  '  Washout  alpha  selection 

ELSEIF  (T  =  2)  THEN 

ALPHA(I)  =  ALPHA2(I):  '  Dead-wgt.  induced  alpha 

'  selection 

ELSEIF  (T  =  3)  THEN 

ALPHA(I)  =  ALPHA3(I):  '  Camber  alpha  section 

ELSEIF  (T  =  4)  THEN 

ALPHA(I)  =  ALPHA4(I):  '  Aileron  float  selection 

ELSEIF  (T  =  5)  THEN 

ALPHA(I)  =  ALPHA(I)  +  ALPHA1 ( I )  +  ALPHA2 ( I )  +  ALPHA3 ( I ) 

+  ALPHA4(I):  '  Complete  wing  selection 
END  IF 

XIN(I)  =  ALPHA(I) 
FOR  J  =  1  TO  KMAX 

A(I,  J)  =  ASYM(I,  J)  -  S(I,  J)  *  Q(I):  NEXT  J 
NEXT  I 

'  Find  Spanload  Soln.  L/Q  from  XIN(KMAX) 

GOSUB  2000 

'  S/R  returns  XIN(I)  as  L/Q  vector  of  length  KMAX 

LIFT  =  0 ! :  MOMENT  =  0  !  :  FOR  I  =  1  TO  KMAX 

LIFT  =  LIFT  +  XIN(I) 

MOMENT  =  MOMENT  +  XIN(I)  *  (XMAC25  -  .5  *  (X1(I)  + 

X2 ( I ) ) ) :  NEXT  I 
'  Find  Lift  &  Moment  Coefficient 
CL  =  LIFT  *  2!  *  DELB  /  S 

CM  =  MOMENT  *  2!  *  DELB  /  (S  *  MAC) 

NP  =  .25  -  (CM  /  CL) 

CLTOT  =  CLREQD  -  CL 

PRINT  "CL  =";  CL;  ■ ,  CM  =»;  CM;  •• ,  CLTOT  =»;  CLTOT. 

'PRINT  "Neut.  Pt.  at  (%  Cmac)";  NPs  PRINT 

GOSUB  1100:  '  Print  Header 

CAVE  =  .5  *  CR  *  (1!  +  TR) 

FOR  I  =  1  TO  M:  '  Sum  up  L/Q  at  eta  value  to  get  CLC 

ETA  =  (I  -  .5)  /  M:  CLC  =  0! 

FOR  J  =  1  TO  N:  K=  (I  -  1)  *N+J 

CLC  =  CLC  +  XIN(K) 

NEXT  J 

'  Find  struct,  twist  due  to  airloads  at  front  panel  for 

'  sta.  "I" 

TWIST  =0!:J2=(I-1)*N+1 

FOR  K  =  1  TO  KMAX 

TWIST  =  TWIST  +  S(J2,  K)  *  Q(I)  *  XIN(K) :  '  Struct  twist 

'  due  to  airload 

NEXT  K 
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C3  =  CR  *  (1!  -  (C0NST1  *(I-.5)/M)):  '  Chord  at 

'  control  point  sta. 

CLSECT  =  CLC  /  C3 

CLC  =  CLC  /  (CL  *  CAVE) 

GOSUB  1210:  '  Print  eta,  cLc/CLCave,  cL(sect)  &  struct. 

'  twist 

NEXT  I 

GOTO  990 

900  '  Anti-Symmetric  Solution  Branch 

'  Set  up  alpha  for  roll  damping  due  to  pB/2V  =  HELIX 
INPUT  "Enter  run  desired:  0=Roll  helix, 

l=Aileron  deflection,  RUN=",  R 
'  Calculate  available  Pb/2V  (rad)  at  specified  velocity 
'  (kts) 

HELIX  =  (.21009  -  .0008744  *  VEL  +  .000001  *  VEL  *  2) 
IF  (R  =  1)  GOTO  905 
FOR  I  =  1  TO  M 
ALPDAMP  =  (I  -  .5)  /  M 

'  Multiply  alpha  for  roll  damping  by  available  roll 
'  helix  angle  (rad)  at  given  velocity  (kts) 
ALPDAMP  =  ALPDAMP  *  HELIX 
FOR  J  =  1  TO  N 

K=(I-1)*N+J:  XIN(K)  =  ALPDAMP:  NEXT  J 
NEXT  I 
GOTO  908 

905  '  Use  following  statement  for  aileron  control 
'  effectiveness 

FOR  I  =  1  TO  KMAX:  XIN(I)  =0!:  NEXT  I 

'  Calculate  aileron  deflection  (rad)  available  at  given 
'  velocity  (kts) 
DELAIL  =  -(90.30699  -  .65754  *  VEL  +  .001771  *  VEL  *    2    - 

.0000016  *  VEL  A  3)  *  .01745 
'  Deduct  float  angle  from  available  control  wheel 
'  aileron  deflection 
FLAIL  =  -(1.4548  -  .0090025  *  VEL  +  .0000444  *  VEL  ~  2) 

*  .01745 
DELAIL  =  DELAIL  -  FLAIL 

XIN(39)  =  .2  *  DELAIL:  XIN(44)  =  .2  *  DELAIL 
XIN(49)  =  .2  *  DELAIL 
XIN(35)  =  .35  *  DELAIL 
XIN(40)  =  DELAIL:  XIN(45)  =  DELAIL:  XIN(50)  =  DELAIL 

908  '  Solve  Spanload  Problem 

FOR  I  =  1  TO  KMAX:  FOR  J  =  1  TO  KMAX 

A(I,  J)  =  AANT(I,  J)  -  S(I,  J)  *  Q(I):  NEXT  J 

NEXT  I 

'  Use  S/R  ELU/SLVB  to  solve  anti-symm  L/Q 

GOSUB  2000 

ROLL  =  0 ! :  FOR  I  =  1  TO  KMAX 
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ROLL  =  ROLL  +  XIN(I)  *  Y3 ( I ) :  NEXT  I 

CROLL  =  -ROLL  *  2!  *  DELB  /  (S  *  B) 

IF  (R  =  0)  THEN 

PRINT  "Roll  Damping  Deriv.  CLP  =" ;  CROLL  /  HELIX 

ELSE 

PRINT  "Aileron  Control  Deriv.,  CLa  =  ",  CROLL  /  DELAIL 

END  IF 

910  'IF  Q  >  0  THEN  GOTO  915 
'CROLLQ0  =  CROLL 

915  '  Dummy  Skip  Statement 

'EROLL  -  CROLL  /  CROLLQ0 

'GOSUB  1200s   '  Print  results 

'  Branch  for  cLc,CLsect, Twist 

GOSUB  1100;  '  Print  Header 

CAVE  =  .5  *  CR  *  (1!  +  TR) 

FOR  I  =  1  TO  M   'Sum  up  L/Q  at  ETA  to  get  CLC 

ETA  =  (I  -  .5)  /  M:  CLC  =  0! 

FOR  J  =  1  TO  N:  K=  (I  -  1)  *N+J 

CLC  =  CLC  +  XIN(K) 

NEXT  J 

'  Find  struct,  twist  due  to  airloads  at  front  panel  for 

'  sta.  "I" 

TWIST  =0!:J2=(I-1)*N+1 

FOR  K  =  1  TO  KMAX 

TWIST  =  TWIST  +  S(J2,  K)  *  Q(I)  *  XIN(K) 

'  Struct  twist  due  to  airload 

NEXT  K 

C3  =  CR  *  (1!  -  (CONST1  *  (I  -  .5)  /  M)):  '  Chord  at 

'  control  point  sta. 

CLSECT  =  CLC  /  C3 

'CLC  =  CLC  /  (CL  *  CAVE) 

GOSUB  1210:  '  Print  eta,  cLc,  cL(sect)  &  struct,  twist 

NEXT  I 

'GOSUB  1100 

'CAVE  =  .5  *  CR  *  (1!  +  TR) 

'FOR  I  =  1  TO  M:   '  Defaulted  for  N  =  1  (Chrdwse  sta.) 

'ETA  =  (I  -  .5)  /  M 

'CLC  =  XIN(I)  /  CAVE:  GOSUB  1200 

'NEXT  I 


990  '  Q  =  Q  +  1! 
'NEXT  MM 
'GOTO  700 


1000  END 
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1100  '  S/R  to  Print  header,  whole  selection  shown 

PRINT  TAB(3);  "Mach  No.1';  TAB (14);  "CLalpha";  TAB (25); 

"CmCl" ;  TAB (35);  "N.P.";  TAB (45);  "Clp" 
PRINT  TAB(5);  "Q,  psi";  TAB(14);  "CLalpha";  TAB(25); 

"CmCl";  TAB(35);  "  N.P.";  TAB(45); 
"Elift";  TAB (55);  "Del-NP" 
PRINT  TAB (10);  "Roll  Control  Evaluation,  Mid  Ail.  at 

0.05" 
PRINT  TAB(5);  "Q,  psi";  TAB(15);  "Croll";  TAB(25); 

"Eroll" 
PRINT 

PRINT  TAB(6);  "2Y/B";  TAB(13);  "Clc/CLcave" ;  TAB(26); 
'  "CI";  TAB (35);  "Twist";  TAB (45);  "Degrees" 

RETURN 
1200  '  S/R  to  Print  results,  whole  selection  shown 

'PRINT  USING  "#####.####";  Q;  CL;  CM;  NP;  ELIFT;  DELNP 
'PRINT  USING  "#####.####";  Q;  CROLL;  EROLL 
1210  PRINT  USING  "#####.####";  ETA;  CLC;  CLSECT;  TWIST; 

TWIST  /  .01745 
WRITE  #1,  ETA,  CLC,  CLSECT,  TWIST,  TWIST  /  .01745 
RETURN 

2000  '  S/R  ELU 

'  Tri-Diagonalizes  the  input  matrix  A(KMAX,KMAX) 
'  Input  Column  vector,  XIN(KMAX),  gets  replaced  by 
'  the  output  which  is  returned  to  the  main  program 

NM1  =  KMAX  -  1 

FOR  K  =  1  TO  NM1:  KP1  =  K  +  1 

FOR  I  =  KP1  TO  KMAX 

G  =  -A(I,  K)  /  A(K,  K):  A(I,  K)  =  G 

FOR  J  =  KP1  TO  KMAX 

A(I,  J)  =  A(I,  J)  +  G  *  A(K,  J):  NEXT  J 

NEXT  I 

NEXT  K 

GOSUB  2050:   '  Use  S/R  SLVB  for  next  step 

RETURN 

2050  '  S/R  SLVB 

'  Solves  the  Tri-Diagonalized  matrix  [A]  obtained 

'  from  ELU  by  back  substitution 

NM1  =  KMAX  -  1:  NP1  =  KMAX  +  1 

FOR  K  =  1  TO  NM1 

KP1  =  K  +  1 

FOR  I  =  KP1  TO  KMAX 

XIN(I)  =  XIN(I)  +  A(I,  K)  *  XIN(K) 

NEXT  I 

NEXT  K 

XIN(KMAX)  =  XIN(KMAX)  /  A(KMAX,  KMAX) 

FOR  K  =  2  TO  KMAX 

I  =  NP1  -  K 

Jl  =  I  +  1 
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FOR  J  =  Jl  TO  KMAX 

XIN(I)  =  XIN(I)  -  A(I,  J)  *  XIN(J) 

NEXT  J 

XIN(I)  =  XIN(I)  /  A(I,  I) 

NEXT  K 

RETURN 

3000  'S/R  to  Determine  [S]  Structural  Infl.  Coeff.  matrix 

3010  '  Generate  [My]  and  [Mx]  matrices,  order  M  x  KMAX 
'  First  get  diagonal  elements 
K  =  1:  FOR  I  =  1  TO  M:  FOR  J  =  1  TO  N 
MY(I,  K)  =  .5  *  (XEA(I)  -  X1(K)  -  .75  *  DELB  * 

TAN(SWP(K) ) ) 
MX(I,  K)  =  .25 
K  =  K  +  1:  NEXT  J 
NEXT  I 

'  Now  get  Upper  triangular  elements 
VALUE  =  1 ! :  FOR  L  =  1  TO  M  -  1 
K  =  L*N+1:K2=M-L 
FOR  I  =  1  TO  K2:  FOR  J  =  1  TO  N 

MY(X,  K)  =  XEA(I)  -  XI (K)  -  .5  *  DELB  *  TAN(SWP(K)) 
MX(I,  K)  =  VALUE:  K  =  K  +  1:  NEXT  J 
NEXT  I 

VALUE  =  VALUE  +  1 ! :  NEXT  L 
'  Test  circuit 

'  PRINT  "Print  out  of  MX(I,J)» 
'  FOR  I  =  1  TO  M:  FOR  J  =  1  TO  KMAX 
'  PRINT  USING  "##.##";  MX(I,  J); 
'  NEXT  J 
'  PRINT  :  NEXT  I 

3020  '  Generate  [DELA(I,J)],  order  KMAX  x  M 
K  =  1 

FOR  I  =  1  TO  M:  FOR  J  =  1  TO  N 
DELA(K,  I)  =  1! 
K  =  K  +  L:  NEXT  J 
NEXT  I 

'  Test  circuit 

'  PRINT  "Print  out  of  DELA(I,J)" 
'  FOR  I  =  1  TO  KMAX:  FOR  J  =  1  TO  M 
'  PRINT  USING  "##.##";  DELA(I,  J); 
'  NEXT  J 
'  PRINT  :  NEXT  I 


3030  '  Generate  [U(I,J)]  sguare  matrix,  order  M  x  M 
FOR  I  =  1  TO  M:  FOR  J  =  1  TO  M 
IF  I  =  J  THEN  U(I,  J)  =  .5 
IF  I  <  J  THEN  U(I,  J)  =  0! 
IF  I  >  J  THEN  U(I,  J)  =  1! 
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NEXT  J 

NEXT  I 

'  Test  circuit 

'  FOR  I  =  1  TO  M:  FOR  J  =  1  TO  M 

'  PRINT  USING  "##.##";  U(I,  J); 

'  NEXT  J 

'  PRINT  :  NEXT  I 


3040  '  Generate  UGJ(I,J)  &  UEI(I,J)  matrices,  order  M  x  M 
'  Embed  DELB~2  type  of  constants 
DELB2  =  DELB  A  2:  TANDELB2  =  TANEA  *  DELB2 


FOR  1=1 
'  DELB2  = 
UGJ(I,  J) 
UEI(I,  J) 
NEXT  I 


TO  M:  FOR  J  =  1  TO  M 

1!:  GJ(J)  =  6!  -  J  :  '  for  test  circuit 

=  DELB2  *  U(I,  J)  /  GJ(J) 


=  TANDELB2  *  U(I,  J)  / 
=  TANDELB2  *  U(I,  J)  / 
'Test  Circuit 

'  PRINT  "Sample  of  UGJ(I/J)" 
'  FOR  I  =  1  TO  M:  FOR  J  =  1  TO  M 
'  PRINT  USING  "##.##";  UGJ(I,  J); 
'  NEXT  J 
'PRINT  :  NEXT  I 


EI(J) 
EI(J) 


NEXT 
NEXT 


J 
J 


3050  '  Start  gathering  terms  for  [S]  matrix  assembly 

DELSIN  =  DELB  *  SINEA:  DELCOS  =  DELB  *  COSEA 

FOR  I  =  1  TO  M:  FOR  J  =  1  TO  KMAX 

F(I,  J)  =  0! :  G(I,  J)  =  0! 

FOR  K  =  1  TO  M 

F(I,  J)  =  F(I,  J)  +  UGJ(I,  K)  *  (COSEA  *  MY(K,  J)  + 
DELSIN  *  MX(K,  J) ) 

G(I,  J)  =  G(I,  J)  +  UEI(I,  K)  *  (SINEA  *  MY(K,  J)  - 
DELCOS  *  MX(K,  J) ) 

NEXT  K 

NEXT  J 

NEXT  I 
3060  '  Form  the  [S]  matrix,  order  KMAX  x  KMAX 

FOR  I  =  1  TO  KMAX:  FOR  J  =  1  TO  KMAX 

S(I,  J)  =  0! 

FOR  K  =  1  TO  M 

S(I,  J)  =  S(I,  J)  +  DELA(I,  K)  *  (F(K,  J)  +  G(K,  J)) 

NEXT  K 

NEXT  J 

NEXT  I 

'  [S]  matrix  is  now  available 

RETURN 
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APPENDIX  B 


C      P-3  LEADING  EDGE  FAILURE  ANALYSIS 
C      TWO  DIMENSIONAL  PANEL  METHOD 

PROGRAM  PANEL 

C THIS  FORTRAN  PROGRAM  ITERATIVELY  COMPUTES  THE  PRESSURE 

C  AND  FORCE  DISTRIBUTION  AROUND  THE  P-3  AIRFOIL  AT  24 

C  SELECTED  WING  STATIONS  AND  WRITES  A  LOAD  FILE  FOR  USE 

C  IN  MSC  PAL2  FINITE  ELEMENT  ANALYSIS  SOFTWARE.   INPUT 

C  CONSISTS  OF  THE  AIRFOIL  COORDINATES  (X,Y)  FROM  A  FILE 

C  CALLED  P3.DAT,  THE  24  WING  STATION  LOCATIONS  IN 

C  INCHES  AND  THE  NODE  NUMBERS  FROM  THE  FINITE  ELEMENT 

C  MODEL  AT  WHICH  FORCES  AND  CONSTRAINTS  WILL  BE  APPLIED. 

C  OUTPUT  CONSISTS  OF  THE  PAL2  LOAD  FILE  (LE.LD),  A 

C  TABULAR  FILE  IDENTIFYING  THE  COEFFICIENT  OF  PRESSURE, 

C  DYNAMIC  PRESSURE  AND  FORCES  AT  A  GIVEN  PANEL  MID  POINT 

C  ON  THE  LEADING  EDGE  OF  THE  AIRFOIL.   SECTION  LIFT 

C  COEFFICIENT  (CI)  IS  MATCHED  TO  AN  INPUT  Cl  REQUIREMENT 

C  (GENERATED  BY  AN  AEROELASTIC  SPAN  LOAD  ANALYSIS 

C  PROGRAM)  THROUGH  ITERATION  OF  ANGLE  OF  ATTACK. 

PARAMETER (N=53,M=N+1 ,PI=3. 14159265 ,RHO=. 0023769, 
+  PAMB=14.7) 

REAL  X(M) ,Y(M) , THETA(M) , BETA(M,M) , R(M,M) ,XM(M) ,YM(M) , 
+      A( 100 ,101) ,AN(M,M),BN(M,M),AT(M,M),BT(M,M), 
+      SUMBN(M,M) ,SUMB1,SUMB2,SUMA(M) ,SUMB(M) ,VTAN(M) , 
+      Q(M) ,CP(M) ,P(M) ,NUM(M) ,DEN(M) , ALPHA, V,VEL,AOA, 
+      CY(M),CX(M) ,CM(M),FY(M),FX(M) , VTOT(N) , CYT, CXT, CMT 
+      FYT,FXT,CD,CL,CLRQD(100) ,WS(100) ,XML(M) ,YML(M) , 
+      LOADM(M) ,INTWIST,OUTWIST,TWIST,XTOPDIS,XBOTDIS, 
+      XDISTOP(48) ,XDISHIN(24) ,ZDISTOP(48) ,ZDISHIN(24) 

INTEGER  K,I, J, NODE ( 100 , 100 ) ,NTOP(100) ,NHIN(100) 

OPEN ( UNIT=5 1 , FILE= ' P3 . DAT ' , STATUS= ' UNKNOWN ' ) 
OPEN ( UNIT=5  7 , FILE= ' DUMP . DAT ' , STATUS= ' UNKNOWN ' ) 
OPEN ( UNIT=58 , FILE= ' WS . DAT ' , STATUS = 'UNKNOWN ' ) 
OPEN ( UNIT=60 , FILE= ' LE . LD ' , STATUS= ' UNKNOWN ' ) 
OPEN ( UNIT=6 1 , FILE= ' TOP . DAT ' , STATUS= ' UNKNOWN ' ) 
OPEN ( UNIT=62 , FILE= ' HINGE . DAT ' , STATUS= ' UNKNOWN ' ) 
OPEN ( UNIT=63 , FILE= ' NODESa . DAT ' , STATUS= ' UNKNOWN ' ) 
OPEN ( UNIT=64 , FILE= ' NODESb . DAT ' , STATUS= ' UNKNOWN ' ) 
OPEN ( UNIT=6  5 , FILE= ' NODESc . DAT ' , STATUS= ' UNKNOWN ' ) 
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C READ  IN  AIRFOIL  DATA 

DO  10  K=1,N 

READ(51,500)  X(K),Y(K) 
10  CONTINUE 
500  FORMAT (2F10. 6) 

X(N+1)=X(1) 
Y(N+1)=Y(1) 

C READ  IN  FINITE  ELEMENT  NODE  NUMBERS  FOR  LATER  USE 

C      (READ  AS  NODE (PANEL  NO ., RIB)) 

DO  12  1=17,36 

READ(63/575)  (NODE ( I , J) , J=l , 8 ) 
READ(64,575)  (NODE( I , J) , J=9,  16 ) 
READ(65,575)  (NODE ( I , J) , J=17 , 24 ) 
12  CONTINUE 
575  FORMAT(9I5) 

C..... WRITE  HEADER  FOR  FINITE  ELEMENT  LOAD  FILE 

WRITE (60, 590) 

C REQUEST  AND  RECEIVE  AIRSPEED  INPUT 

PRINT*, 'ENTER  VELOCITY  IN  KNOTS  (XXX. XX) ' 

READ*,VEL 

V  =  VEL*1.6878 

C REQUEST  AND  RECEIVE  FRONT  SPAR  TWIST  ANGLES 

PRINT*, 'ENTER  TWIST  AT  2Y/B=.35  ( rad) ( .XXXX) ' 

READ*,TWIST35 

PRINT*, 'ENTER  TWIST  AT  2Y/B=.45  (rad) ( .XXXX) ' 

READ*,TWIST45 

PRINT*, 'ENTER  TWIST  AT  2Y/B=.55  ( rad) ( .XXXX) ' 

READ*,TWIST55 


C CALCULATE  CHORD  LENGTH  OF  P-3  WING  AT  SELECTED  WING 

C      STATIONS 


C SET  UP  OUTER  LOOP  TO  ITERATIVELY  CALCULATE  PRESSURES 

C      FOR  A  SERIES  OF  WING  STATIONS 


99 


C INPUT  WING  STATIONS,  Cl  AND  LOAD  MULTIPLIER  FROM  FILE 

C      WS.DAT  AND  CALCULATE  CHORD  LENGTH  AT  THE  WING  STATION 

502  FORMAT(F7.1,F8.4,F6.2) 
DO  15  L=l,24 

READ (58,502)  WS(L), CLRQD ( L )  , LOADM ( L ) 
PRINT* , ' ITERATION= ' , L 
PRINT*, 'WS=',WS(L) 
PRINT* , ' CLRQD= ' , CLRQD ( L ) 

CHORD  =  227.*(1-(0.6*(2.*WS(L)/1188. ) ) ) 
PRINT* , ' CHORD= ' , CHORD 

C CALCULATE  THE  THICKNESS  FRACTION  AT  A  GIVEN  WING 

C      STATION  (USE  TO  LINEARLY  VARY  AIRFOIL  THICKNESS  WITH 
C      WING  STATION) 

DTHICK  =  1-((1-(12./14. ))*(2.*WS(L)/1188. )) 
PRINT*, 'THICKNESS  RATIO  =', DTHICK 

C SCALE  THE  THICKNES  OF  THE  AIRFOIL  AT  EACH  WING  STATION 

C      BY  THE  THICKNESS  FRACTION 

DO  16  K=1,N 

Y(K)=Y(K)*DTHICK 
16  CONTINUE 

C SET  STARTING  ALPHA 

ALPHA  =  (PI/180. ) 
C. SET  COUNTER  START  FOR  AOA  ITERATION 

KOUNT  =  1 

C. BEGIN  ITERATION  OF  ALPHA  TO  MATCH  CLREQD  =  CL 

1100  CONTINUE 

C FIND  NORMALIZED  MID-POINT  OF  PANEL 

DO  30  1=1, N-l 

XM(I)=0,5*(X(I)+X(I+1) ) 
YM(I)=0.5*(Y(I)+Y(I+1) ) 

C FIND  MID-POINT  OF  PANEL  IN  INCHES 

XML ( I ) =XM ( I ) *CHORD 
YML ( I ) =YM ( I ) *CHORD 
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C CALCULATE  THETA  (ANGLE  BETWEEN  PANEL  AND  CHORD  LINE) 

NUM(I)=Y(I+1)-Y(I) 
DEN(I)=X(I+1)-X(I) 
THETA(I)=ATAN2(NUM(I) ,DEN(I) ) 

C CALCULATE  R  (DISTANCE  BETWEEN  PANEL  MID-POINTS) 

DO  20  J=1,N 

R(I,J)=SQRT( (XM(I)-X(J) )**2+(YM(I)-Y(J) )**2) 

C CALCULATE  BETA   (ANGLE  AT  MID  POINT  OF  PANEL  I, 

C      ENCLOSING  END  POINTS  OF  PANEL  J) 

BETA(I,J)=ATAN2( ( ( YM( I )-Y( J+l ) ) * (XM( I )-X( J) ) 

#  -((XM(I)-X(J+1)))*(YM(I)-Y(J))), 

#  (<XM(I)-X(J+1))*(XM(I)-X(J)) 

#  +(YM(I)-Y(J+1))*(YM(I)-Y(J)))) 


20     CONTINUE 
30  CONTINUE 

C. CALCULATE  AUGMENTED  "A"  MATRIX 

SUMB1=0.0 
SUMB2=0.0 
DO  50  1=1, N-l 

SUMBN ( I , N ) =0 . 0 
DO  40  J=1,N-1 

IF  (I.EQ.J)  THEN 
AN(I,J)=0.5 
AT (I, J) =0.0 
BN(I, J) =0.0 
BT(I,J)=0.5 
ELSE 

AN(I, J)=(SIN(THETA(I)-THETA(J) ) *LOG(R( I , J+l )/R( I , J) ) 

#  +COS(THETA(I)-THETA( J) )*BETA(I/ J) )/(2*PI) 

BN ( I , J )  =  ( COS ( THETA ( I ) -THETA ( J ) ) *  LOG (R(I,J+1)/R(I,J)) 

#  -SIN(THETA(I)-THETA( J) )*BETA(I, J) )/(2*PI) 
AT(I,J)=-BN(I,J) 

BT ( I , J ) =AN ( I , J ) 
END  IF 

SUMBN(I/N)=SUMBN(I,N)  +  BN(I,J) 
IF  (I.EQ.l)  THEN 

SUMB1=SUMB1  +  BT(1,J) 
ELSEIF  (I.EQ.N-1)  THEN 

SUMB2=SUMB2  +  BT(N-1,J) 
END  IF 
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A(I,J)  =  AN(I,J) 

A(I,N)  =  SUMBN(I,N) 

A(N,J)  =  AT(1,J)  +  AT(N-1,J) 

A(N,N)  =  SUMB1  +  SUMB2 
40     CONTINUE 
50  CONTINUE 

DO  60  1=1, N-l 

A(I,N+1)  =  -V*SIN(ALPHA-THETA(I) ) 
60  CONTINUE 

A ( N , N+ 1 ) =-V* ( COS ( ALPHA-THETA ( 1 ) ) + 
+  COS ( ALPHA-THETA ( N- 1 ) ) ) 

CALL  GAUSS (N, A) 

GAMMA=0 . 0 
DO  80  1=1, N 
Q(I)=A(I,N+1) 
GAMMA=A(N,N+1) 
80  CONTINUE 

DO  100  1=1, N-l 

SUMA(I)=0.0 

SUMB(I)=0.0 
DO  90  J=1,N-1 

SUMA( I ) =SUMA( I ) +Q ( J ) *AT ( I , J ) 

SUMB ( I ) =SUMB ( I ) +BT ( I , J ) 

VTOT ( I ) =SUMA ( I ) +GAMMA*  SUMB ( I ) +V*COS ( ALPHA-THETA ( I ) ) 

VTAN ( I ) =VTOT ( I ) / V 

CP(I)=1.0-(VTAN(I) )**2.0 

P(I)=( (CP(I)*0.5*RHO*(V**2. ) )/144. ) 
90  CONTINUE 
100  CONTINUE 

C CALCULATE  FORCE  COEFFICIENTS  (CX,CY)  AND  FORCES  (FX= 

C      CHORDWISE,  FY=NORMAL  TO  CHORD)  AT  PANEL  MID  POINT 

CYT=0 
CXT=0 
CMT=0 
FXT=0 
FYT=0 
DO  110  1=1, N-l 

CY(I)=-1.0*CP(I)*(X(I+1)-X(I) ) 

CX(I)=CP(I)*(Y(I+1)-Y(I) ) 

C MULTIPLY  BY  CHORD  FOR  FULL  SCALE  SOLUTION 

FY( I )=CHORD*CY( I ) *0 . 5*RHO* ( V**2 . )/144 
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FX(I)=CHORD*CX(I)*0.5*RHO*(V**2. )/144 
CM(I)=CP(I)*( (X(I+1)-X(I) )*XM(I)+ 
+         (Y(I+1)-Y(I) )*YM(I) ) 
CYT=CYT+CY(I) 
CXT=CXT+CX(I) 
FYT=FYT+FY(I) 
FXT=FXT+FX(I) 
110  CONTINUE 

C GET  FORCE  TOTALS  FORWARD  OF  FRONT  SPAR  (.15C) 

FYNT=0 . 0 
FXNT=0 . 0 
DO  115  1=18,35 

FYNT=FYNT+FY(I) 

FXNT=FXNT+FX(I) 
115  CONTINUE 

C CALCULATE  APPROXIMATE  TOTAL  FORCES  ON  COMPLETE  LE 

C      SECTION  (BETWEEN  NACELLES)  BY  TAKING  FORCES  AT 
C      MID  POINT  *  92  INCHES 

IF  (L.EQ.7)  THEN 

FYTOT=FYNT*92 

FXTOT=FXNT*92 
END  IF 

C .....  PERFORM  AXIS  ROTATION  TO  FIND  CL 

CL=CYT* COS ( ALPHA )-CXT* SIN (ALPHA) 

C PERFORM  ITERATION  OF  ANGLE  OF  ATTACK 

DIFF=.0001 

CLDIFF=CLRQD ( L ) -CL 

KILL=50 

DELTA=0.1745 

IF  (KOUNT.GT.KILL)  THEN 

PRINT* , ' COUNTER  EXCEEDED ' 

GOTO  111 
ELSEIF  (CLDIFF.GT.DIFF)  THEN 

PRINT'  (I3,2F10.6)  ',KOUNT,  ALPHA ,  CL 

KOUNT=KOUNT+l 

ALPHA=ALPHA+ ( DELTA*ABS ( CLDIFF ) ) 

GOTO  1100 
ELSEIF  ( CLDIFF. LT. 0.0)  THEN 

PRINT' (I3,2F10.6) ' ,KOUNT, ALPHA, CL 

KOUNT=KOUNT+l 

ALPHA=ALPHA-( DELTA* ABS( CLDIFF) ) 

GOTO  1100 
ELSE 
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PRINT' (4F10.5) ' , ALPHA, CL,FXNT,FYNT 

WRITE(57,570)  VEL, WS ( L) , ALPHA, CLRQD ( L) , CL , FXNT , FYNT 
WRITE(57,560) 
END  IF 

C... WRITE  PRESSURES  AND  FORCES  TO  OUTPUT  FILES 


DO  120  1=17,36 
WRITE (57, 580) 
+ 
120  CONTINUE 


I,XM(I),XML(I),YML(I),CP(I),P(I), 
FX(I)*LOADM(L) ,FY( I ) *LOADM(L) 


C. . . . .WRITE  LOADING 
APPLICATION ( PAL2 ) . 


FILE  (LE.LD)  FOR  FINITE  ELEMENT 
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FX(17)-FX( 
FY(17)=FY( 
FX(23)=FX( 
FY(23)=FY( 
FX(25)=FX( 
FY(25)=FY( 
FX(26)=FX( 
FY(26)=FY( 
FX(28)=FX( 
FY(28)=FY( 
FX(30)=FX( 
FY(30)=FY( 
FX(36)=FX( 
FY(36)=FY( 
DO  125  1=1 

J=I 
IF  (I.LE.2 

WRITE (60 

ELSEIF  (I. 
WRITE (60 

ELSEIF  (I. 
WRITE (60 

ELSEIF  (I. 
WRITE (60 

ELSEIF  (I. 
WRITE (60 

END  IF 
CONTINUE 


17)*0c5 

17)*0.5 

23)+0.5*FX 

23)+0.5*FY 

25)+0.5*FX 

25)+0.5*FY 

26)+FX(27) 

26)+FY(27) 

28)+0.5*FX 

28)+0.5*FY 

30)+0c5*FX 

30)+0.5*FY 

36)*0.5 

36)*0.5 

7,36 


(24) 
(24) 
(24) 
(24) 


(29) 
(29) 
(29) 
(29) 


3)  THEN 

,591)  NODE(J,L 

NODE(J,L 
EQ.25)  THEN 
,591)  NODE(J,L 

NODE ( J , L 
EQ.26)  THEN 
,591)  NODE(J,L 

NODE ( J , L 
EQ.28)  THEN 
,591)  NODE(J,L 

NODE ( J , L 
GE.30)  THEN 
,591)  NODE(J,L 

NODE(J,L 


,FX(I)*LOADM(L), 
,FY(I)*LOADM(L) 

,FX(I)*LOADM(L), 
,FY(I)*LOADM(L) 

,FX(I)*LOADM(L) , 
,FY(I)*LOADM(L) 

,FX(I)*LOADM(L), 
,FY(I)*LOADM(L) 

,FX(I)*LOADM(L) , 
,FY(I)*LOADM(L) 
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560  FORMAT (2X, 'PANEL' , 3X, 'XM( I ) ' , 4X, 'XML(I) ',4X, 'YML(I) ', 
+        5X, 'CP(I) ',5X, 'P(I) ',6X, 'FX(I) ' ,5X, 'FY(I) ' ) 

570  FORMAT(//,2X, 'VEL= ' , F10 . 1 , / , 2X, 'WS= ' , F10 . 3 , / , 2X, 

+      'ALPHA=',F10.5,/,2X, 'CLRQD= ' , F10 . 5 , /, 2X, 'CL=' , 
+      F10.5,/,2X,  'FXNT=',F10.5,/,2X,  'FYNT=  '  ,  F10  .  5  ,  /  ) 

580  FORMAT(2X,I3,7F10.5) 

590  FORMAT( IX, 'FORCES  AND  MOMENTS  APPLIED  0') 

591  FORMAT (IX, 'FX ' , 17 , F10 . 5 , /, IX, 'FZ ' , 17 ,F10 . 5 ) 
111  CONTINUE 

15  CONTINUE 

C CALCULATE  TWIST  DISPLACEMENTS  OF  FINITE  ELEMENT  NODES 

INTWIST=0.82155*(TWIST45-TWIST35)+TWIST35 
OUTWIST=0. 89371* (TWIST55-TWIST45 )+TWIST45 
TWIST=OUTWIST-INTWIST 
XTOPDIS=TWIST*8 . 8 
XBOTDIS= (-1.0) *XTOPDIS 

DO  130  1=1,24 

J=I-1 

XDISTOP(2*I-l)=(XTOPDIS/23. )*J 

XDISTOP(2*I)=(XTOPDIS/23. )*J 

XDISHIN(I)=(XBOTDIS/23)*J 

TT=(TWIST/23. )*J 

XT=(XTOPDIS/23. ) * J 
IF  (TT.LE.0625)  THEN 

ZDISTOP( 2*1 ) =0 . 5*TT*XT 
ELSE 

ZDISTOP (2*I)=(-0.5)*( TT- . 125 ) *XT 
END  IF 

ZDISTOP(2*I-1)=(-0.5)*TT*(XT+1.0) 

ZDISHIN( I )=0 . 5*TT*XT 
130  CONTINUE 

WRITE (5 7, 5 55)  FXTOT,FYTOT 
555  FORMAT (/,2X, 'APPROX  TOTAL  FX( 92" )= ' ,F10 . 1, 
+        /,2X, 'APPROX  TOTAL  FY( 92" )= ' ,F10 . 1 ) 

C READ  AND  WRITE  NODES  AND  CORRESPONDING  BOUNDARY 

C      CONDITIONS  IN  FINITE  ELEMENT  LOAD 

WRITE (60, 597) 

READ(61,592)  (NTOP( I ) , 1=1, 48 ) 

READ(62,601)  (NHIN( I ) , 1=1, 24 ) 

WRITE(60,593)  (NTOP( I ) , 1=1, 48 ) 

WRITE (60, 5931)  (NTOP( I ) ,XDISTOP( I ) ,1=1,48) 
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C ENFORCE  Y-AXIS  DISPLACEMENT  CONSTRAINT  ALONG  FRONT 

C      EDGE  OF  SPAR  CAP  FLANGE  ONLY.   (RELEASE  BACK  EDGE  TO 
C      SIMULATE  SINGLE  LINE  OF  SCREW  FASTENERS.) 

WRITE (60, 5932)  (NTOP(2*I) ,1=1,24) 

WRITE (60 ,5933)  (NTOP(I) , ZDISTOP( I ) , 1=1, 48 ) 

WRITE(60,594)  (NHIN( I ) ,XDISHIN( I ) ,1=1,24) 

WRITE(60,599)  (NHIN( I ) , 1=1, 24) 

WRITE (60 ,600)  (NHIN( I ) , ZDISHIN( I ) , 1=1, 24) 

WRITE(60,595)  (NHIN( I ) , 1=1, 24 ) 

WRITE(60,596)  (NHIN( I ) , 1=1, 24 ) 

WRITE (60, 598) 

592  FORMAT(I5) 
601  FORMAT(I5) 

593  FORMAT(lX, 'RA',I7,1X, '0' ) 

5931  FORMAT(lX, 'TX',I7,1X,F8.4) 

5932  FORMAT(lX, 'TY',I7,1X, '0' ) 

5933  FORMAT(lX,  'TZ  ',I7,1X,F8.  4) 

594  FORMAT(lX, 'TX' , 17 , 1X,F8 . 4) 

599  FORMAT(lX, 'TY',I7,1X, '0') 

600  FORMAT(lX, 'TZ ' , 17, 1X,F8 . 4) 

595  FORMAT(lX, 'RX',I7,1X, '0' ) 

596  FORMAT(lX, 'RZ',I7,1X, '0') 

597  FORMAT(lX, '--  BLANK  LINE  — ',/, 

+        IX, 'DISPLACEMENTS  APPLIED  0') 

598  FORMAT(lX,' —  BLANK  LINE  — ',/, 
+        IX, 'SOLVE',/, IX, 'QUIT') 

STOP 
END 


SUBROUTINE  GAUSS (N, A) 

C THIS  SUBROUTINE  PERFORMS  GAUSS  ELIMINATION  WITH 

C      PIVOTING. 

INTEGER  PV  !  PIVOT  INDEX. 

DIMENSION  A( 100, 101) 
EPS=1.0 
10  IF  (1.0+EPS.GT.1.0)  THEN 
EPS=EPS/2.0 
GOTO  10 
END  IF 
EPS=EPS*2.0 
EPS2=EPS*2.0 
DO  1010  1=1, N-l 

PV=I 
DO  J=I+1,N 
IF  (ABS(A(PV,I) ) .LT.ABS(A(J,I) ) )  PV=J 
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END  DO 

IF  (PV.EQ.I)  GOTO  1050 

DO  JC=1,N+1 

TM=A(I, JC) 

A(I, JC)=A(PV, JC) 

A(PV, JC)=TM 
END  DO 
1050  IF  (A(I/I) .EQ.O)  GOTO  1200   !  KICK  OUT  IF  SINGULAR 
DO  JR=I+1,N         !  ELIMINATION  OF  BELOW  DIAGONAL. 
IF  (A(JR,I).NE.0.0)  THEN 

R=A(JR,I)/A(I,I) 
DO  KC=I+1,N+1 

TEMP=A(JR,KC) 

A(JR,KC)=A(JR,KC)  -  R*A(I,KC) 
IF  (ABS(A(JR,KC) ) .LT.EPS2*TEMP)  A( JR,KC)=0 . 0 

C IF  THE  RESULT  OF  SUBTRACTION  IS  SMALLER  THAN  2  TIMES 

C      EPS  TIMES  THE  ORIGINAL  VALUE,  IT  IS  SET  TO  ZERO. 

END  DO 

END  IF 
1060  END  DO 
1010  CONTINUE 

IF  (A(N,N).EQ.O)  GOTO  1200 

A(N,N+1)=A(N,N+1)/A(N,N) 

DO  NV=N-1,1,-1      !  BEGIN  BACKWARD  SUBSTITUTION. 

VA=A(NV,N+1) 
DO  K=NV+1,N 

VA=VA-A(NV,K)*A(K,N+1) 
END  DO 

A(NV,N+1)=VA/A(NV,NV) 
END  DO 
RETURN 

1200  PRINT*, 'MATRIX  IS  SINGULAR' 

STOP 
END 
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